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algorithms for analyzing the non-stationary data are 
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turbulent boundary layer include the non-stationary mean 
velocity, turbulence intensity, power spectral density and the 
autocorrelation and cross correlation functions. Since the 
data are periodic, both ensemble averaging and special non- 
stationary empirical model techniques are employed. Each of 
the algorithms has been implemented in the FORTRAN language 
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l. INTRODUCTION 


The need to understand the nature of turbulent fluid flow 
has long been a field of aerodynamic research. Much of the 
productive work in the field has been experimental in nature. 
Since turbulence is a stochastic or random process, analysis 
of experimental data requires the use of statistical methods. 
Due to the complex nature of turbulent fluid flow and its most 
frequent applications, most of the experimental work has been 
restricted to steady (statistically stationary) turbulence. 
Statistical methods for stationary random processes are well 
understood. Recent investigations by Holmes, Obara and Yip 
[Ref. 1] and Howard [Ref. 2] have begun to extend experimental 
research into unsteady (non-stationary) turbulent fluid flow 
conditions. The most recent work by Howard consists of an 
investigation of the time varying characteristics of the 
boundary layer of an airfoil immersed in ae propeller 
slipstream. These experiments require the application of non- 
stationary statistical data analysis methods. 

Future wind tunnel tests have been proposed by Howard 
[Ref. 2] to continue investigations into the time dependent 
nature of the boundary layer of an airfoil subjected to 


periodic turbulent wake disturbances. The proposed 


experimental work has motivated the need for non-stationary 
statistical data analysis methods with a sound mathematical 
and statistical basis. 

This thesis describes an effort to develop a complete data 
analysis system to support the proposed wind tunnel tests. 
The data Dee may also be applied to other unsteady 
turbulent flows, such as internal flows in turbomachinery. 
Specific algorithms for analyzing non-stationary turbulent 
boundary layer velocity data are identified, developed and 
implemented. Where alternate algorithms are identified, each 
algorithm is evaluated and the best method is recommended 
based on specified criteria. 

Chapter II presents a brief overview of turbulent fluid 
flow and boundary layer theory. Chapter II includes a 
description of the specific statistical quantities used to 
characterize turbulence, a brief description of the results 
of the recent tests by Howard [Ref. 2] and the detailed 
objectives and design constraints for the unsteady turbulence 
data analysis system. Chapter III briefly reviews the basic 
mathematical and statistical concepts which were applied to 
develop the data analysis methods. Chapter IV presents the 
detailed explanations of each algorithm and includes test case 
examples used to evaluate the operation and suitability of the 
computer programs. Conclusions and recommendations are 


presented in Chapter V. Appendix A presents some derivations 


and analyses not included in Chapter IV. Appendix B contains 
the user's guide and FORTRAN source code for each of the 


computer programs. 


Il. NATURE OF THE PROBLEM 


The purpose of this chapter is to present background 
material on the topic of turbulent fluid flow and boundary 
layer theory, to briefly review recent work completed on the 
behavior of unsteady boundary layers in a_e propeller 
slipstream, and then to state the specific objectives of this 
thesis. The following references were found helpful in the 
preparation of this chapter: Bradshaw [Ref. 3], Cebeci and 
Smith [Ref. 4], Frost and Moulden [Ref. 5] and Howard [Ref. 


Ae 


A. WHAT IS TURBULENCE? 


The phenomenon of fluid turbulence occurs in nearly all 
fluid flow situations. To paraphrase the definitions of 
turbulence presented in Cebeci and Smith [Ref. 4}: 


Turbulent fluid motion is an irregular fluid flow 
condition which, in general, appears when fluids flow past 
solid surfaces or when neighboring streams of fluid pass 
over one another in shear layers. Turbulence is 
characterized by a random high frequency variation in 
various quantities (e.g., velocity, pressure, temperature, 
etc.) with time and space coordinates, so that 
statistically distinct average values can be discerned. 


Turbulent fluid flow is a "mixing" flow where the fluid 
streamlines are intersecting, twisting and curling about each 


other. In contrast, laminar flow is characterized by smooth, 


parallel streamlines with no mixing and no random fluctuations 
in the flow parameters. 

Bradshaw [Ref. 3] defines the study of turbulence as the 
art of understanding the Navier-Stokes equations without 
actually solving them. The Navier-Stokes equations, which 
have been derived from Newton's second law of motion, are the 
governing equations of fluid flow and are three dimensional 
non-linear partial differential equations in four independent 
Variables (three space dimensions and time). These equations 
express complicated relationships between the fluid quantities 
and have proven resistant to explicit solution, resulting in 
"too many unknowns and not enough equations". The primary 
deterrence to the solution of the Navier-Stokes equations has 
been the lack of a suitable mathematical turbulence model. 

Therefore, much of the productive study of turbulent fluid 
flow results from statistical analysis of empirical (measured) 
data to characterize the nature of turbulence. The velocity 
(or other parameters of interest) of a turbulent fluid in 
motion is often modeled by defining the instantaneous velocity 
aS a sum of a mean velocity u and a randomly fluctuating 


velocity, u', with a zero mean, as in the equation: 


u(t) = u + u’ (t) (II-1) 
Instantaneous Mean Fluctuating 
Velocity Velocity Velocity 


where u(t) is the velocity in the direction of the mean flow 


(the x direction). The y and z turbulent velocity components 


are similarly defined using the variables v and w respectively 
except that the mean velocities v and w are zero. In the 
context of statistics, steady turbulence is defined as a fluid 
flow with a constant mean velocity and a constant variance 
high frequency fluctuating velocity component of constant 
variance (statistically stationary). Unsteady turbulence is 
characterized as a fluid flow with time-dependent mean 
velocity and/or time-dependent variance fluctuating velocity 
component (non-stationary). Figure II-1 [from Bradshaw, Ref. 
3] presents an illustration of steady vs. unsteady turbulent 
flows. In contrast to turbulent fluid flow, laminar flow 
(either steady or unsteady) is defined as a flow in which the 
instantaneous velocity 1S essentially equal to the mean 
velocity (i.e., the fluctuating velocity component is zero or 


very small: on the order of 0.02% of the mean velocity). 


B. BOUNDARY LAYERS 


A boundary layer is a thin region of fluid flow next to 
a solid surface such as an airfoil in which strong viscous 
effects exist. Due to fluid viscosity, the layer of fluid 
particles touching the surface are attached to the surface and 
have no velocity relative to the surface. The fluid velocity 
increases with distance from the surface until the velocity 
equals that outside the boundary layer. The way in which the 


velocity increases across the boundary layer is called the 
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Figure II-1. EXAMPLES OF STEADY AND UNSTEADY TURBULENT FLOWS 
(BRADSHAW, REF. 3]. 


boundary layer velocity profile (see Figure II-2). The 
boundary layer typically grows in thickness as the fluid flows 
across the surface from the leading edge. 

Boundary layer velocity profiles have been divided into 
different categories based on their shape. The two most 
important profiles are called laminar and turbulent velocity 
profiles and each profile has its own characteristic shape. 
Figure II-3 presents examples of typical laminar and turbulent 
boundary layer velocity profiles. Due to the mixing and 


momentum transfer in a turbulent boundary layer, the turbulent 
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Figure II-2. TRANSITION FROM A LAMINAR BOUNDARY LAYER TO A 
TURBULENT BOUNDARY LAYER ON A SMOOTH FLAT PLATE. 


boundary layer profile is characterized by a higher velocity 
near the surface compared to a laminar profile. The higher 
velocity near the surface causes a higher frictional drag on 
the surface compared to the laminar boundary layer. Boundary 
layers near the leading edge of smooth surfaces usually start 
out laminar. Viscous forces cause boundary layers to 
transition from laminar to turbulent at some distance back 
from the leading edge. One of the methods of reducing 
aerodynamic drag is to delay the transition of the boundary 
layer from laminar to turbulent for as long as possible on the 
entire surface. The boundary layer can be sensitive to 
surface roughness or to external disturbances in the flow, 
such as turbulent wakes from control surfaces or propellers. 
Such disturbances can cause early transition of the boundary 


layer from the laminar to the turbulent case. 
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Figure II-3. LAMINAR VERSUS TURBULENT BOUNDARY LAYER PROFILE. 
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C. PHYSICAL QUANTITIES OF INTEREST 


This section presents some of the physical quantities used 
to describe the characteristics of a turbulent flow. These 
statistical parameters are defined in this chapter using the 
terminology of the aerodynamicist such as in Cebeci and Smith 
[Ref. 4]. Chapter III then reviews the statistics and 
mathematics of random processes which describe many of the 
Same parameters in the terminology of the statistician. 

1. Turbulence Intensity 

Turbulence intensity 1S a measure of the relative 
energy of the turbulence. Turbulence intensity is defined as 
the root mean square of the fluctuating velocity component 


where the mean square value is given by: 
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The mean square value of the v and w velocity components are 
defined similarly. If the fluid flow is steady, then this 
quantity (which is similar to the statistical variance) is 
constant and is computed using time averages. If the flow is 
unsteady, then the turbulence intensity varies with time and 
must be computed from ensemble averages (see Chapter III). 
Turbulence intensity is normalized by dividing by the mean of 
the velocity at the outer edge of the boundary layer or some 
other suitable mean velocity. For the x, y and z velocity 


components, turbulence intensity are: 
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The magnitude of the normalized turbulence intensity parameter 
is on the order of 0.01 (1.0%) for highly turbulent flows and 
on the order of 0.0001 (0.01%) for laminar flows [Ref. 4}. 
The turbulence intensity of the free stream flow in a typical 
wind tunnel is 0.3%. In specially designed low turbulence 
wind tunnels, turbulence intensities of 0.02% are attained. 
2. Turbulence Length Scale 

A measure of the characteristic length scale of the 
turbulent flow (i.e., the average size of a turbulent eddy) 
is obtained by measuring the spatial correlation of the 
velocity between two points in the flow. This requires two 
or more velocity sensors at different locations in the flow. 
If the fluctuating velocity component, u , is small compared 
to the mean velocity, then the eddies do not’ change 
Significantly as they pass a given point. This common 
Situation allows the use of time correlation (autocorrelation) 
to characterize the turbulence length scale instead of spatial 
correlation (eliminating the need for multiple sensors). 
Autocorrelation is defined using the  aerodynamicist 
terminology as: 


mage) u’ (t+7) = Lim a S (ute) eulten (Gt 7, ) (Ti-4) 


lend 
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where for steady turbulence i is a time index and for unsteady 
turbulence i 1S an ensemble index. For three dimensional 
measurements, autocorrelations of v and w may also be 
defined in a Similar manner. 

Though harder to understand physically, cross 
correlations are also important to the study of turbulent 
flows. These are a measure of the correlation between two 

/ y 


different components of fluctuating velocity, u’ and v’° for 


example. The cross correlation is defined by: 


Sas oes ee N 

u(t) v/ (tt) = lim q uy (t)v," (t47) (II-5) 
ee eee et N 

vi (t) uf (t+r) = lim gy S v,(t)u’,(t+7) 


N7o@ 


it 
— 


where for steady turbulence 1 is a time index and for unsteady 
turbulence i is an ensemble index. The cross correlation of 
u andv (or any two components) for zero lag are Known as 
the "Reynolds stress" terms, and are important to 
understanding the convective interaction between the 
turbulence and the mean flow [Ref. 4]. 
3. Frequency Spectra 

With respect to turbulence, frequency spectra are 
measures of the power contained in the turbulent fluctuating 
velocity components at a given frequency. The various 


frequency spectra (auto-spectra and cross-spectra) contain the 


2 


same information contained in the various’ correlation 
measures, but transformed into the frequency domain. 


Frequency spectra are defined in more detail in Chapter III. 


D. PREVIOUS UNSTEADY TURBULENCE INVESTIGATIONS 


The boundary layer on a surface such as an airfoil can, 
in many situations, be sensitive to external disturbances in 
the flow. One example of such an external disturbance has 
been investigated experimentally by Howard [Ref. 2]. This 
investigation involves the study of the boundary layer of an 
airfoil immersed in a propeller slipstream (see Figure II-4). 
The propeller wake subjects a given point on the airfoil to 
periodic turbulent pulses. The propeller wake provides a 
convenient method to study the effect on a boundary layer of 
an external turbulent disturbance. The periodic nature of the 
propeller wake pulses allow the rapid collection of many 
repeated disturbances for ensemble averaging, which is 
necessary due to the unsteady nature of the fluid flow. 

"Single-wire" hot wire anemometer velocity data were 
collected and boundary layer velocity profiles were 
constructed using a simple ensemble average of time averaged 
segments (see Chapter IV). A single-wire velocity sensor 
measures the magnitude of the velocity normal to the wire 
(that is qu2 +- - Howard [Ref. 2] concluded that the wing 


boundary layer was found to vary from the laminar to turbulent 
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Figure II-4. AN AIRFOIL IN A PROPELLER SLIPSTREAM [Howard, 
Ref. 2]. 


condition in a periodic manner due to the passage of the 
propeller wake at a frequency of 20 to 40 Hz [see Figure II- 
5 from Howard, Ret. 2i- The length of time the turbulence 
remains in the boundary layer is a strong function of the 
pressure gradient (the change in pressure with location on the 
surface in the direction of mean flow). The turbulence within 
the boundary layer persists longer as the pressure gradient 
varies from favorable (dP/dx < 0) to adverse (dP/dx>0). Te 
was also observed that the shape of the velocity profile has 
a strong dependence on the instantaneous level of external 


turbulence. Furthermore, the instantaneous drag coefficient 
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of the airfoil section was found to vary from that expected 
for a laminar boundary layer to below that expected for a 
conventional fully turbulent boundary layer. In other words, 
the mean drag was well below that expected for a fully 
turbulent boundary layer. No correlation or spectral analyses 


were performed. 
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Figure II-5. EXAMPLES OF PERIODIC BOUNDARY LAYER OSCILLATIONS 
(Howard, Ref. 2]. 

An understanding of the time-dependent behavior of an 
airfoil boundary layer responding to external disturbances 


will provide basic information to studies of the use of 


dynamic maneuvering for increased fighter aircraft agility. 
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Dynamic maneuvering involves the use of large deflections of 
aircraft control surfaces at very high angles of attack (post- 
stall angles of attack) in order to achieve maximum agility 
particularly at low speeds. Dynamic maneuvering creates the 
possibility of an aircraft wing or control surface being 
exposed to the turbulent wake disturbance from other surfaces. 
Another potential application is to the extremely turbulent 
flow (on the order of 20%) in turbomachinery such as jet 


engines. 


E. THESIS OBJECTIVES 


1. Future Unsteady Turbulence Investigations 

Future wind tunnel tests have been proposed for fall 
1988 to continue the study by Howard [Ref. 2]. The 
investigation will be expanded to include "dual-wire" hot wire 
anemometer measurements which will allow resolution of the u 
and v velocity components within the boundary layer. Boundary 
layer profiles and turbulence intensity measurements will be 
determined for the two component velocity data. Also, 
autocorrelation analyses (turbulence scales), cross 
correlation analyses (Reynolds stresses) and spectral analyses 
will be accomplished using both single and dual wire velocity 
sensors to examine these properties in a periodically 


disturbed boundary layer. Of interest are the differences 
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between a conventional turbulent boundary layer and one 
excited by freestream turbulence. 
2. Thesis Objectives 

This thesis describes an effort to develop a set of 
data analysis algorithms and software (with appropriate 
mathematical justification) to support the unsteady turbulence 
investigation described in the section above. The work 
includes investigations into alternate methods for determining 
some of the statistical properties used to characterize the 
unsteady (non-stationary) periodically perturbed boundary 
layer. Where alternate algorithms are developed, the 
strengths and limitations of each method are examined and the 
"best" algorithm is recommended based on criteria such as 
accuracy, error sensitivity, speed of execution and ease of 
implementation. 

Specifically, data analysis algorithms are required 
to extract certain statistical parameters used to characterize 
the nature of turbulent fluid flow from the "raw" velocity 
data. The specific unsteady turbulence parameters required 
to support the planned wind tunnel tests described above are 
listed as follows: 

Required Turbulence Parameters 
* Mean Velocity 
* Fluctuating velocity component 


* Turbulence intensity 


7 


- Spectral (frequency) characteristics 
¢ Auto/cross correlation (turbulence scales) 
The basic design requirements for the completed data 
analysis software include the following constraints: 
Software Design Constraints 
¢ Software must run on an IBM PC/AT class microcomputer. 
- Software must be written in FORTRAN. 


¢- Data processing will be performed after testing (no real 
time processing requirements). 


¢ Input file formats are defined by the analog-to-digital 
conversion system. 
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lll. MATHEMATICAL BACKGROUND AND 
DEVELOPMENT 


The purpose of this chapter is to present a brief overview 
of the basic mathematical concepts which were applied to 
develop the specific data analysis methods for the unsteady 
turbulence investigation. The important concepts include the 
statistics of random processes associated with time series 
analysis such as time and ensemble mean and variance, 
autocorrelation and time cross correlation, probability 
distributions and stationary vs. non-stationary random 
processes. Also important is the Discrete Fourier Transform 
(DFT), particularly as it applies to spectral estimation. The 
following references were used extensively in the preparation 
of this chapter; Bendat and Piersol [Ref. 6], Brigham [Ref. 
7), Strum and Kirk [Ref. 8], Hamming [Ref. 9], Marple [Ref. 


10], and Otnes and Enochson [Ref. 11]. 


A. STATISTICS OF RANDOM PROCESSES 


1. Definition of Terms 
Turbulent fluid motion is characterized by random 
fluctuations of the observed quantities (such as velocity or 


pressure) with respect to space and time. Random data are not 
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deterministic and as such, cannot be described by explicit 
mathematical expressions. Statistical analysis methods are 
used to extract and summarize useful information that 
characterizes the nature of the random process. A sequence 
of observations or measurements is called a "sample function", 
a "time series" or a "time history record". Physically, 
turbulent motion is a continuous’ process. But the 
measurements sample the continuous process at discrete times. 
In addition to being discretized in time, the observed 
quantity (e.g. velocity) is also quantized for storage in a 
digital computer. For example, given a 12 bit analog to 
digital conversion, the measured velocity may only take on 
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or 4096 distinct values within some predefined range. All 
of the time series considered here will be discrete and 
quantized and are denoted as: 

U, OF U(t,), I=H eye ae 
A collection of sample functions is called an "ensemble" and 
is denoted as: 

u, (t; ) 

An ensemble of all possible sample functions defines the 
underlying “random process". Figure III-1 illustrates an 
ensemble of N time history records. 
2. Probability Distribution Functions 


Given a discrete random variable, u,, there exists a 


finite probability that u, will be equal to a particular value 


Zo 





Figure III-1. AN ENSEMBLE OF TIME HISTORY RECORDS [BENDAT AND 
PIERSOL, REF. 6]. 


at time index i. For example, if u is the fluid velocity at 


a point in the turbulent boundary layer and the index i 
represents a particular instant in time (t,), then there is a 
probability value (between zero and one), say 0.3 for this 
example, that u will be equal to a particular velocity, U, 
where U is one of 4096 distinct values allowed for u by its 
digital representation. This is expressed as: 
p(u,) = PROB(u,=U) = 0.3 (uit. 7) 

The probability mass function, p(u,;) of a discrete random 
variable defines how the total probability of 1.0 is 


distributed over the range of possible values of the random 
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variable. The probability mass function is usually not known 
for an empirically measured set of random data. 

Given a continuous random variable, u(t), there exists 
a finite probability that u(t) will be between two particular 
values U and U +AU at a given time t. For the continuous 
case, u may take on any value. The probability density 
function, f(u), is defined by the differential relation as 


QU approaches zero: 


PROB (U<u(t) <U+AU) 


Lar Tli-2 
f(u) = lim Xt (Ii 
Note that, 
J. £(u) du=1 (III-3) 


For many practical applications, a continuous random process 
which has been quantized for a digital computer (such as the 
velocity example above) will still be discussed in terms of 
its continuous probability density function even though it 
actually takes on only discrete values. This is a valid 
approximation as long as the quantization is sufficiently 
"fine". (At least 512 to 1024 distinct values over the range 
of the measured quantity will be considered sufficient.) 
There are many theoretical probability distribution 
functions or "distributions" which have been defined for both 
discrete and continuous random variables. The most common is 


the continuous "Gaussian" or "Normal" distribution. If a 
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random variable u has mean #, and variance oes (see definitions 
below), then it is normally distributed if its probability 


density function is given by: 





(u-p)? 
p(u) = —— Exp | -——+— (III-4) 
270, 26, 








The normal distribution appears often in nature and 
much of its importance is due to the Central Limit Theoren, 
which states that the distribution of a sufficiently large 
linear combination of finite variance independent random 
variables approaches a normal distribution. This holds even 
if the distribution of each member of the independent linear 
combination is not normally distributed. This is important 
because an empirically measured random process can often be 
modeled as a linear combination of many independent random 
variables (say 20 or 30 as a rule of thumb). Then by the 
Central Limit Theorem, the distribution of the measured random 
process may be assumed normal. If a random variable is 
normal, then the mean and variance provide a complete 
description of the nature of the random process. Bendat and 
Piersol [{Ref. 6] presents a good summary of many other 
continuous and discrete probability density functions and 
their uses, as well as a more detailed discussion of the 


Central Limit Theoren. 
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3. Mean Value 
The most basic statistic is the mean value. The mean 
value is also known as the expected value or average value. 
In terms of the discrete random variable u and the probability 
mass function p,(u), the mean value of the ensemble of time 
history records, u,(t;) at a fixed time t, is defined by the 


equation: 


y (ti) = ELC) ] =Ly,(t,) p(w) (III-5) 


where the "E[ ]" notation is the expected value of u. The 
expected value operator may also be applied to any function 


of a random variable and is defined by: 


E[g(u,)] = x g(u,) p(u) (III-6) 


In terms of the probability density function f(u), and the 
continuous random variable u(t), the mean (expected) value of 
the ensemble of time history records at fixed time t is 


defined by: 


E[u,(t)] = fut) f(u) du (Itt=7) 


For practical measurements, the number of sample functions 
(and the length of time) are finite. The sample mean (an 
unbiased estimate of the true mean) of a finite ensemble of 
N time history records at a fixed time t, 1s definedeby: 


A = Lae 
H(t) = B(t)) = yz y(t) (III-8) 


u 
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where the notation denotes an estimator of the true 
parameter. If a random process is a member of the special 
class of “ergodic" random processes (see below), then the 
ensemble averaged mean value of the random process is 
invariant with time and is equal to the time averaged mean 


value of a single (infinite) time history. The time mean of 


a single time series is defined by: 


yO) = Lim # Ee w(t) = w(t) (III-9) 
u +0 ie 


For a practical finite time history record, the time sample 
mean (again, an unbiased estimate of the true mean for ergodic 


random data) is given by the equation: 


s u, (t,) (III-10) 


BOO = Gk) = fe > 


4. Variance 
Another important statistic is the "variance" of the 
random process. The variance 1s a measure of the mean square 
deviation or dispersion of the data about its mean. The 
square root of the variance is called the standard deviation. 
The variance of an ensemble of discrete time history records 


is defined by the equation: 


2 
g(t.) = E fu, (t,) -H, (t,) | | (III-11) 





= F(acty-v(t)}) pc] 


k= 
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where the mean value, u,, is computed by equation III-5 and 
o, is called the standard deviation. For a practical fame. 
ensemble, an unbiased estimate of the true variance is given 


by the sample variance: 


aA? 1 
6. (Ol) =_S (Ce a 


IMI 


2 
uw, (t,) -4(t,) | (Mii) 


pr 
au 


Similar to the discussion of mean above, if a random process 
is ergodic (see below) then the ensemble variance is equal to 
the variance of a single infinite time history record and is 


constant with time. The time variance is defined by: 
; 1 ©™ 2 , 
O° (k) = Lim yay 2 (Sh) | crate) (III-13) 


where the mean value is computed by equation III-9. For a 
practical finite time history record, an unbiased estimate of 


the time variance is given by the time sample variance: 


A 2 2 1 M — .2 
6, (k) = S2(k) = ge ¥ (u,(t,)-i,) (III-14) 
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5. Autocorrelation and Cross Correlation Functions 
Correlation is a measure of the degree of linear 
relationship that exists between two random variables, u and 
v, at a given instant in time. The correlation function 
introduces another variable, 7, which is a time lag between 
the two random variables. The correlation function 
essentially gives a correlation measure for each of the 


possible values of the lag *% (that is, the amount of 
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correlation when one time series leads or lags the other). 
When u equals v, the correlation function is’ called 
autocorrelation. When u is different than v, the correlation 
function is called cross correlation. 

In terms of the expected value operator, the 


autocorrelation of a real valued time series is defined as: 
R,,(t,,t,+7) = E}u,(t,)u, (t,+7) (III-15) 


For an ergodic random process, the ensemble definition of 
autocorrelation is equal to the time averaged definition for 
a single infinite duration time series. The time averaged 
definition of autocorrelation is: 


N 
Pet) = Lim a UES Re sere) (III-16) 
1 


= © =-N 


An unbiased estimate of the autocorrelation function of a 
real, discrete, finite, ergodic random process with zero mean 


u(t,) 1S given by: 

A ak N-r 
Ri (7) = yop 2 ul(t,)u(tyt7) xr=0,1,2,...,m (III-17) 

j=] 


where r, the lag number, is defined by T such that T= rat, 
At is the sample interval and m is the maximum lag number. 
The maximum lag number, m, can go up to N-1, but for typical 
computations is usually much less than N (on the order of 1/10 
N). This equation is expressed only for positive lag numbers. 


However, it has been shown that the autocorrelation function 
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is an even function; thus R, (7) equals R,(-%). See Bendat 
and Piersol [Ref. 6]. 
For two random processes u and v, the cross 
correlation function is defined by: 
R,,(t,,t,+7) = E[u, (t,) v, (t,+7) | (III-18) 
Ryy (ty, ty+7) =E| vy (t,) uy (t,+7) | (Tit 
An unbiased estimate of the cross correlation function for a 


practical finite data set consisting of two ergodic random 


variables with zero mean is given by: 


N=r 
A -20 
Rw(7) = fer > UCen(eie (Ti teees 
T=rAt 
r=0., U2 je 
- N-r 
Ru(7) = oe Sv (ti)u(tit7) (III-21) 
T=ypAt 
P=0 2G ee 


The cross correlation function estimate is in general neither 
odd or even, but for negative lag numbers, it can be shown to 
satisfy the relation that, for real valued data, eG equals 
R,(-f) The autocorrelation function is a special case of the 
cross correlation function for u equal v. The auto 
correlation and cross correlation functions may be estimated 
directly using equations III-17 and III-20, III-21 or may be 
computed indirectly using Fourier transform methods which are 
computationally more efficient for large data sets. Fourier 


transform methods are addressed in section III-B. 


28 


6. Spectral Density Functions 
The autospectral density function is a measure of the 

energy content of a time history record at a given frequency. 
The classical definition of the autospectral density function 
of a discrete random variable u, is the Fourier transform of 
the autocorrelation function. In terms of the discrete 
Fourier transform, the autospectral density function is 
defined by: 

Syy(£,) =Ot, S Ryy (it) EXP[-J27iK/N] (11-22) 


k=0,11,72,... 


where f, is a particular frequency, t is the sample interval 
and I is the sample index such that =i tt. By the Wiener- 
Khinchine relations [Ref. 6] it has been proved that the 
corresponding spectral density functions and correlation 
functions constitute Fourier transform pairs. For a real 
valued time series, equation III-22 may be simplified to yield 


the "single sided" spectral density function: 
Guy (f,) =40t D Ryy (iMt) cos (27ik/N) (III-23) 
oes... 
which is only defined for positive frequencies. Equation III- 
23 1s therefore twice the magnitude of the even function "two 
Sided" spectral density function given by equation III-22. 
For practical finite length time series, substitute the 
autocorrelation function estimate (equation III-16) into 


equation III-23. The importance of this topic warrants a more 
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detailed discussion. Additional material will be included in 
the review of the Fourier transform methods in Section III-B. 
7. Classes of Random Data 

Random processes may be classified as either 
stationary or non-stationary (see Figure III-2). Stationary 
Gata are further classified as ergodic or non-ergodic. Non- 
stationary data are classified as non-stationary mean, non- 
stationary variance, non-stationary frequency and other 
specialized classes. A random process is non-stationary (the 
most general case) when the ensemble averaged mean and 
autocorrelation function (and other ensemble properties) vary 
with time. The special classes of non-stationary data, non- 
stationary mean, variance and frequency, imply that the mean, 
variance or frequency content of the ensemble vary with time, 
but the other ensemble averaged properties are time invariant. 
In the most general case of non-stationary data, "everything" 


varies with time. 


| STATIONARY WON- STATIONARY | 


— I 
NON - ERGODIC 


ERGODIC WON-STATIONARY 
L_________! MEAN 


NON -STATIONARY NON - STATIONARY 
FREQUENCY OTHER 





Figure III-2. RANDOM DATA CLASSIFICATIONS [BENDAT AND PIERSOL, 
REF. 6]. 
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A random process is classified as "weakly stationary" 
when the ensemble mean and the autocorrelation function are 
time invariant. "Strongly stationary" refers to the special 
class of data where all ensemble averaged statistical 
properties (moments) are time invariant. 

A random process is classified as “ergodic" if it is 
stationary and if the time averaged mean value _ and 
autocorrelation function (and all other time averaged 
properties) are equal to the corresponding ensemble averaged 
values. Thus for an ergodic random process, it is possible 
to derive desired statistical information about the entire 
random process from a single time history record of sufficient 
duration. The mathematics associated with the analysis of 


random processes is much simpler if the process is ergodic. 


B. DISCRETE FOURIER TRANSFORM METHODS 


1. Definition of Terms 

A physical process (in general, either a random or 
deterministic process represented by real or complex numbers) 
is typically described in the "time domain". For example, the 
turbulent velocity measurements u(t) where the velocity is 
described as a function of time (i.e., the typical time 
history record). A physical process may also be described in 
the "frequency domain" where the velocity is given by its 


amplitude, U(f), as a function of its frequency, f. In 
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general, the amplitude is also a complex number. U(f) and 
u(t) are essentially two different representations of the same 
data record. Both representations contain exactly the same 
information regarding the process, presented differently. A 
method for transforming between these two representations is 
the "fourier transform" and the “inverse fourier transform". 

Fourier transforms are defined for both continuous 
and discrete processes. The planned wind tunnel tests 
described in Chapter II involve only discretized data. 
Therefore, the discussion here is limited to discrete 
(sampled) data. Moreover, most modern time series analyses 
are accomplished on discretized data to take advantage of the 
power and speed of digital computers. Some specific issues 
associated with sampled data will be discussed before 
reviewing the discrete fourier transform. 

2. Data Sampling 

For a continuous process u(t) which is discretized or 
sampled at equal time intervals, the sampling interval, At, 
1s the amount of time between samples. The sampling rate, 
f., 1s the inverse of the sampling interval and is usually 


$s 


expressed in units of Hertz (Hz) or samples per second (SPS): 


f= a (TET=omy 


The "Sampling Theorem" (see Strum and Kirk [Ref. 8] 


and Brigham [Ref. 7]) states that if a continuous signal, 
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u(t), contains no frequency components above a frequency Lael 
then the signal can be uniquely represented by equally spaced 
samples u(t.) only if the sampling frequency f, 1s greater than 
2 f,,,> Furthermore, the original signal u(t) can be recovered 


exactly from the sampled signal by the formula: 


(III-25) 


= sin] 2mf-(ti-idt) 


where f, is the Nyquist critical frequency defined by: 


s 


2At ~ 2 


Kh 
Oo 
il 
rh 


(III-26) 


Unfortunately, many realizable processes may contain 
a very wide band of frequency components. This may cause the 
required minimum sample rate of 2f,,. to be impractically high. 
The phenomenon of "aliasing" is a result of sampling a 
wideband process. Data aliasing causes frequency information 
at frequencies above the Nyquist frequency to be "folded" 
spuriously down into the frequency range between -f, and f,. 


The Nyquist critical frequency, f,, 1s the folding frequency. 


c! 


In other words, the sample rate, f must be greater than 


g / 
twice the frequency of the highest frequency present in the 
data to preclude aliasing. An alternate method to prevent 
aliasing is to low pass filter the continuous (analog) signal 


before sampling to remove all high frequency components of the 


data above the frequency range of interest. Then, the data 
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are sampled at a frequency of greater than twice the maximum 
frequency of interest. 
3. The Discrete Fourier Transform 
To transform an N point discrete time series Uo jinog 
the time domain to the frequency domain U, (in general, both 
are complex), the discrete Fourier transform (DFT) is defined 


by: 


N 
U, = & u, EXP(-j2nik/N) (III-27) 


ome 
| aad 


Here k is the frequency index running from -N/2 to N/2 (or 0 
to N-1) and N is the total number of points in the time series 
u,. The DFT is said to map u, into U,. If u, is consideredmms 
be a sequence of samples of a continuous function u(t), then 


the DFT times the sampling interval is an approximation of the 


continuous Fourier transformation of u(t). That is, 


N 
U(f,) = At U, = At 2, EXP (-j2mik/N) (III-28) 


where f,=k/N t and f, varies between -f. and f,. U(f,) is an 
approximation of the continuous Fourier transform because the 
Sampled time series is of finite duration and finite frequency 
bandwidth (a sampled and truncated approximation of the 
infinite duration continuous signal). In the development of 
the DFT as a special case of the continuous Fourier transforn, 
it is assumed that the finite time series {uu} 15 a4 Sameae 


period of an infinite real periodic function. This makes the 
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"infinite periodic function" an even function so that by the 
properties of symmetry, the DFT is also an even function. 
That is, U(f) equals U(-f). Also of importance is another 
symmetry property which states that if u, is real and even 
(and if it is a finite length sampled version of a continuous 
function, it must be even), then the DFT will also be real and 
even (see Press, et al, [Ref. 12] and Brigham, [Ref. 7]). 
The inverse discrete Fourier transform (IDFT) is 


defined by the following equation: 


N 
u=% 2U, EXP(j2ik/N) (11 me) 


The IDFT transforms a frequency domain function back into the 
time domain. Given U,, by equation III-29, the IDFT will 
recover exactly the u,'s that were input into the DFT 
equation. Note that the summation is over the frequency index 
k for the IDFT instead of the time index i. Equations III-27 
and III-29 form a Fourier transform pair. Given U(f,) 


(equation III-28) the IDFT is defined by: 
ie te: 
Ui = War 20 (£,) EXP(j271ik/N) (III-30) 


The only difference is the factor of 1/At which causes the 
Fourier transform pair III-20 and III-25 to depend on the time 
scale while the pair III-19 and III-21 do not. The key point 


is to be consistent in the Fourier transform pair used. 
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4. Correlation Functions 
Correlation functions (auto and cross) were defined 
in section III.A.5. above. In this section, the "indirect" 
DFT method for computing the auto or cross’ correlation 
functions is described. The Discrete Correlation Theorem 
Ref. 7] states that the "circular" correlation functions may 


be determined in the frequency domain through the following 


relationship: 
N 
n° _& +7) = + SUV .EXP(j20rk/N) (III-31) 
Stare) aes SS u)) N 2, i Yk 
T=rAt 


PH02) 2 peewee 

In words, the circular correlation function of two discrete 
time series u(t;) and v(t,) may be computed by multiplying the 
complex conjugate of the DFT of the first time series by the 
DFT of the second time series, and then computing the IDFT of 
the product. This is essentially a linear convolution. 
Stated another way, the expression UV, forms a Fourier 
transform pair with the left hand side of equation III-31. 
(Again, when u equals v then equation III-31 yields the 
autocorrelation.) 

The nature of the circular correlation function 
results in the correlation values at positive and negative 
lags being summed. To obtain the standard correlation 
functions (equations III-17, III-20, and III-21) instead of 


the circular correlation, add N zeros to the end of u and v 
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(or at least m zeros, where m equals the maximum lag number 
r). Then compute the circular correlation in the normal way. 
The added zeros will cause the negative lags to be separated 
from the positive lags [{Ref. 6]. This comes at a cost since 
the two time series are now twice as long. The standard 
correlation function estimate is given by: 

R (7) = eR uy (7) m=rAt, r=0,1,2,.-.-,N-1l (III-32) 
This method, while more complicated, is more computationally 
efficient for large data sets because it can take advantage 
of an efficient Fast Fourier transform (FFT) algorithm [Ref. 
6). 

5. Autospectral Density Functions 
A theorem important to the development of Fourier 
transform spectral analysis techniques is Parseval's Theoren. 
Parseval's Theorem states that the energy contained in the 
time domain of a time series u,; 1s equal to the energy 
contained in the signal's frequency domain [Ref. 7}. In terms 
of the discrete transform, Parseval's Theorem is stated as 


follows: 


hk 


N N = 
i= = 


= 


This makes intuitive sense because the information contained 
in the time domain representation of a signal is equivalent 


to that contained in the frequency domain. 
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As described in section III.A.6. above, the classical 
definition of the spectral density function is the Fourier 
transform of the autocorrelation function. An alternative 
method applies the Fourier transform directly to the time 
series to obtain a spectral estimate. When an FFT algorithm 
is used to compute the Fourier transform, the spectral 
estimate is obtained many times faster than the classical 
method. (The FFT algorithm will be discussed in the next 
section.) For an infinite length time series, the spectral 


density function is defined by: 


N Z 
4c S u,EXP(-j27ik/N)| | (III-34) 
{=-N 


e 1 
Suu (Fy) =1am 2 (2NT1) ae 





or 


Syu (£,) =Lim e| cowaryaes (9 (80 


and U(f,) is defined by equation III-28. Bendat and Piersol 
(Ref. 6] have shown that equation III-34 is equivalent to 
equation III-22. For a practical finite length time series, 
the basic "raw" sample spectral density function may be 


estimated using the equation: 


S,,(&) = gael (f) Up (III-35) 


This is the original periodogram power spectral density (PSD) 
estimator. Marple [Ref. 10] shows that this PSD estimator 


will yield unstable PSD estimates which possess a variance of 
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the same order as the power estimate. Several averaging 
techniques have been developed to improve the _ spectral 
estimator. The specific technique employed for this thesis 
will be discussed in the next chapter where details of the 
specific algorithm are presented. 
6. Windows in Spectral Estimation 

Another concept of importance to spectral estimation 
is that of window functions. A typical window function is a 
function which starts at zero and smoothly rises to a peak 
value of one, then smoothly falls back to zero. The need for 
window functions or "windowing" arises from estimating 
autospectral density functions using finite length time 
series. When a Fourier series is used to represent a 
function, an exact representation of the function can be 
obtained, in general, only when an infinite number of terms 
in the Fourier series are retained. When the Fourier series 
is truncated, the function is approximated. Similarly, the 
discrete autocorrelation function values can be thought of as 
the coefficients of a Fourier series representing the spectral 
density function. When a time series (and its autocorrelation 
function) is of finite length (truncation), then the resulting 
spectral density function iS an approximation of the true 
spectral density function due to the truncated Fourier series. 
Window functions are used to improve the approximation of the 


spectral density function. 
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Recall that the spectral density function (equation 
III-34) is defined for time samples running in the limit from 
minus infinity to plus infinity. When a practical finite 
length time series is used (equation III-35), this is 
analogous to the coefficients of an infinite Fourier series 
being multiplied by a function which is zero valued outside 
of the length of the finite sample record (i.e., is zero 
before t = 0 and after t = NAt) and is one during the length 
of the time history record. Such a window function is known 
as the "rectangular window". This usually introduces a 
discontinuity at the beginning and end of the time series as 
the DFT "sees" it. Discontinuous functions are especially 
difficult for a truncated Fourier series to approximate. The 
Fourier coefficients corresponding to the higher frequencies 
are needed for the Fourier series to properly represent the 
discontinuous function. But the high frequency coefficients 
are discarded when the Fourier series is truncated. 

Windowing is designed to improve the spectral density 
function approximation. This 1s accomplished by multiplying 
the input time series by a window function which starts at 
zero, smoothly rises to a peak of one, then smoothly falls 
back to zero. This "smoothes out" the discontinuities at the 
beginning and end of the time series record. Then the high 
frequency coefficients of the Fourier series are smaller in 


magnitude and have a much less significant effect on the 
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Fourier series representation of the spectral density 
function. Thus when the Fourier series is truncated, less 
information is discarded and a better approximation is the 
result. 

Figure III-3 presents an example of the resulting 
spectral estimate of a 3750 Hz sine wave (amplitude of one) 
sampled at 15000 Hz. Thus the 3750 Hz signal is centered on 
a frequency bin (i.e., the time series contains an integer 
number of sine wave periods). The resulting spectral estimate 
is exact because the spectral density function is represented 
by a Fourier series with a single non-zero coefficient. 
Figure III-4 presents the more general case where the time 
history record is not an integer number of sine wave cycles, 
in this case a sine wave of 3700 Hz. The 3700 Hz signal falls 
between two frequency bins. Note that the power contained in 
the sine wave has now "leaked" into nearby frequency bins. 
In fact, all of the frequency bins contain some of the power 
that should appear only at 3700 Hz. 

Figure III-5 presents the spectral estimate of the 
3700 Hz sampled sine wave time series which has_ been 
previously multiplied by a "von Hann" window (equation III- 
G6). Note that leakage has been reduced to a few adjacent 
frequency bins yielding a better approximation to the exact 
spectral density function. The magnitude of the power at each 


frequency has been corrected to account for the effect of the 
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Figure III-3. SPECTRAL ESTIMATE OF 3750 HZ SINE WAVE SAMPLED 


AT 15,000 HZ, NO WINDOW FUNCTION. 
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Figure III-4. SPECTRAL ESTIMATE OF 3700 HZ SINE WAVE SAMPLED 
AT 15,000 HZ, NO WINDOW FUNCTION. 
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Figure III-5. SPECTRAL ESTIMATE OF 3700 HZ SINE WAVE SAMPLED 
AT 15,000 HZ, VON HANN WINDOW FUNCTION. 
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window function. The von Hann window function is given by: 


w(t,)=0.5-0.5cos| Sm] = 1 2p feroieesate) 


where N is the total number of points in the time series. 
There are many other window functions which have different 
frequency response characteristics and are used to optimize 
the spectral estimator for specific purposes. Some of these 
will be discussed in Chapter IV with the discussion of the 
specific algorithms. A paper by Harris [Ref. 13] provides an 
excellent description of the use and characteristics of many 
window functions. 
7. The Fast Fourier Transform 

The fast Fourier transform (FFT) 1s a computationally 
efficient algorithm (actually several algorithms) fon 
computing the discrete Fourier transform given by equation 
III-27. If equation III-27 is implemented explicitly on a 
computer, it requires on the order of N° operations 
(multiplies and adds) to compute the DFT of an N point time 
series. The FFT algorithm arrives at the same result in N 
log,N / 2 operations. For example, a time series of length 
1024 would require over one million operations uSing the 
explicit method. But only about five thousand operations are 
required for the FFT, a factor of about 1/200 fewer 
operations! In simplified terms, the FFT algorithm 


accomplishes this by considering the DFT as a system of 
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equations written in matrix form. When the proper matrix 
factorization is applied, zeros are introduced into the 
factors which allow a reduction in the number of required 
computations. The larger the matrix system, the greater the 
reduction in required operations. See Brigham (Ref. 7] for 


a detailed description of the basic FFT algorithn. 
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lV. DATA ANALYSIS ALGORITHMS FOR UNSTEADY 
TURBULENCE 


Ensemble averaging techniques may be used to determine the 
properties of non-stationary data when many repeated time 
history records are available. A general method does not 
exist for analyzing the properties of all types of non- 
stationary data from single time history records [Ref. 6]. 
For single time history records, special techniques must be 
developed that apply to limited classes of non-stationary 
data. 

The problem addressed by this thesis falls somewhere 
between the two extremes of a single data record vs. many time 
history records. The periodic nature of the data allows 
ensemble averaging techniques to be employed. However, 
practical considerations such as computer memory and data 
storage limitations preclude collecting "many" time history 
records. That is, for the current applications, collecting 
50 time history records is feasible but collecting 500 time 
history records is not. However, 50 time history records may 
not be sufficient to obtain "good" estimates of the 
statistical properties using strictly ensemble averaging 


techniques. 
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This chapter describes the algorithms which were developed 
to estimate the statistical properties of unsteady (non- 
stationary) turbulent fluid flow. Both ensemble averaging and 
special (non-stationary empirical model) techniques have been 
employed, often in combination. Each of the algorithms has 
been implemented in the FORTRAN language, and examples are 


presented using a representative test case data set. 


A. TEST CASE DATA--UNSTEADY TURBULENT FLUID FLOW 


Figures IV-1 through IV-3 present three sample records 
(instantaneous velocity vs. time) from an ensemble of 41 
Sample records recorded at a single location above the surface 
in the boundary layer of an airfoil. Recall that a single 
sample record represents one propeller revolution (called a 
"pbulse") and two propeller blade wake passages. This typical 
ensemble of 41 propeller pulses will be used as the test case 
for most of the unsteady turbulence algorithms presented in 
this thesis. The data were obtained from Howard [Ref. 2] and 
are representative of the type of data which will be collected 
during future tests. These data were measured using a single- 
wire velocity sensor so that the velocity measurements 
represent the magnitude of the velocity vector in a plane 


normal to the wire: 


qu2(t)+v2(t) 
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Figure IV-1. PULSE 1 FROM TEST CASE 1 (AN EXAMPLE ENSEMBLE OF 
41 PROPELLER PULSES). SAMPLE RATE = 15,000 HZ. 
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Figure IV-2. PULSE 2 FROM TEST CASE 1 (AN EXAMPLE ENSEMBLE OF 
41 PROPELLER PULSES). SAMPLE RATE = 15,000 HZ. 
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Figure IV-3. PULSE 41 FROM TEST CASE 1 (AN EXAMPLE ENSEMBLE 
OF 41 PROPELLER PULSES). SAMPLE RATE = 15,000 HZ. 
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The velocity data were sampled at 15,000 Hz. The two-blade 
propeller was rotating at approximately 1200 rpm such that a 
turbulent disturbance (propeller wake) occurs at about 40 Hz. 
For simplicity, the velocity is presented in digital counts 
rather than physical units. The data are plotted with 
straight line segments between each of the data points. The 
free stream (aircraft) velocity was 65 miles per hour. Time 
is presented as sample number with sample period t = 6.666 


Sear? 


seconds. Each pulse contains 696 time samples in this 
example. No two-component data (individual measurements of 
u and v) were acquired during the tests of Reference 2. When 
needed for software check-out, two-component data will be 
simulated. 

The standard input file format for single-wire velocity 
data consists of a sequential file containing alternating 
values of velocity and propeller trigger, starting with 
velocity. The propeller trigger is a signal generated with 
a magnetic pickup from one of the propeller blades. Each 
rotation of the blade generates a sharp voltage spike which 
is recorded and used to identify the beginning of each 


propeller pulse. Time is not recorded explicitly but is 


implied from the sample rate. 
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B. MEAN VELOCITY DETERMINATION 


1. Overview 

Equation II-1 describes a common empirical model of 
fluid velocity, u(t), in a turbulent fluid flow as the sum of 
the mean velocity, u, and the (zero mean) stationary random 
fluctuating velocity component, u'. In the boundary layer 
of an airfoil immersed in a propeller slipstream, the non- 
stationary mean velocity varies with time as does the high 
frequency fluctuating component (i.e. the instantaneous 
velocity is non-stationary in mean). The total instantaneous 
velocity is actually measured. So the problem is to extract 
the time varying mean velocity from the measured instantaneous 
velocity. The classical method for determining the mean of 
non-stationary data is ensemble averaging. Ensemble averaging 
works well if a sufficient number of sample records are 
obtained. For this specific problem, the mean velocity varies 
periodically with the propeller rotation at a known constant 
low frequency relative to the turbulent fluctuations. fThis 
"low frequency non-stationary mean" empirical model of the 
velocity data provides an alternate method of analyzing the 
non-stationary data through low pass filtering [Ref. 6]. 
Figures IV-1 through IV-3 illustrate the "low frequency mean" 
nature of the data and justify the use of a low frequency mean 


empirical model. 
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Several different algorithms to determine the time 
varying mean velocity will be examined in this section. The 
algorithms are based on the concept of ensemble averaging, the 
low frequency mean empirical model, or a combination of the 
two. 

2. Ensemble Averaging Algorithm 


a. Description of the Algorithm 


The ensemble averaging algorithm is a straight 
forward implementation of equation III-8. The ensemble 
average is computed using the following steps: 

¢ Identify the input and output data file names. 


- Separate the input time history file into an ensemble of 
individual pulses. An array (ENSMBL) is used to store 
up to 50 propeller pulses of up to 2048 samples each. 
Read the input data file until the first propeller 
trigger is found, indicating the beginning of a pulse. 
Moving down the first column, read each velocity sample 
into the array. When the next propeller trigger is 
encountered, begin at the top of the next column and 
continue reading the velocity samples into the array. 
Repeat the process until all pulses have been put into 
the array. Figure IV-4 illustrates how the velocity data 
are stored into the ensemble array. In the event that 
each propeller pulse does not contain the same number of 
data points, the columns are truncated to the length of 
the shortest column. 


- A row in array ENSMBL represents data at the same fixed 
time point from all of the pulses as illustrated in 
Figure III-1. Per equation III-3, compute the ensemble 
average by first summing the velocity values across each 
row, then dividing each sum by the number of pulses. 
Store the average of each row in the column array AVG. 
AVG represents the ensemble average of all the pulses at 
each time sample. 


- Write the ensemble average array AVG to the output file. 
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Figure IV-4. FORMAT OF ARRAY ENSEMBLE. 
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b. Implementation and Examples 


The above algorithm was implemented in FORTRAN as 
program ENSEMBL. The source code and user's guide for program 
ENSEMBL is presented in Appendix B. Program ENSEMBL was run 
with the test case input file described in section IV.A. 
Figure IV-5 presents a plot of the ensemble average of the 
test case data. Note that the high frequency turbulent 
fluctuations have been greatly reduced but not eliminated 
compared to Figures IV-1 through IV-3. Qualitatively, the 
remaining high frequency fluctuations in the ensemble average 
may be an indication that 41 pulses is not a sufficient number 
to obtain a reasonable estimate of the mean uSing ensemble 
averaging. This observation provides the motivation for 
investigation into alternate methods of determining the non- 
stationary mean velocity. Computer memory and data storage 
limitations preclude any substantial increase in the volume 
of data, thus alternate methods are required to obtain a 
better estimate of the mean velocity. 

The part of the algorithm which separates the 
input file into an ensemble of individual pulses may introduce 
a small timing error when the length of each pulse is 
truncated to the length of the shortest pulse. For test case 
1, an average of three to four samples (out of 700 samples per 
pulse) were truncated from the end of each pulse. This timing 


error is less than one per cent and is considered negligable. 
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Figure IV-5. ENSEMBLE AVERAGE OF TEST CASE 1. 
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The error is primarily due to the inability to control the 
external test conditions, such as airspeed and propeller RPM 
during the test run. The data of test case 1 were collected 
in an airborne environment where controlling the test 
conditions is difficult. Future tests will be conducted in 
a wind tunnel which will allow the test conditions to be 
controlled more precisely. 
3. Moving Average Smoothing Algorithm 


a. Description of the Algorithm 


The moving average smoothing algorithm represents 
a variation on the method used by Howard [Ref. 2] to reduce 
the unsteady turbulence data presented in Reference 2. This 
data analysis technique involves applying simple moving 
average smoothing to the time history data prior to performing 
ensemble averaging. A moving average iS a time average over 
a short segment or "window" of the time history record. The 
time average of the short segment becomes the "smoothed" value 
of the velocity at the center point of the segment. Thus the 
length of the segment,L, is usually specified as an odd number 
of points. As the moving average window is moved across the 
time history record(of length N, where N is much greater than 
L), a smoothed velocity value is computed for each new center 
point of the window. If the moving average window is L points 


long, then the first smoothed value which can be computed is 
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the (L/2 + 1)'th point. The first L/2 points (and the last 
L/2 points) in the time history record are set to zero. 

The assumption associated with the moving average 
is that the time history data are stationary over the short 
window. An empirical model was proposed above for this 
unsteady turbulence data which supposes that the mean is time 
varying at a low frequency relative to the random fluctuating 
velocity component (i.e., the high frequency turbulence). As 
long as this empirical model is a reasonable description of 
the actual data, then the assumption of stationarity over 
short time segments should be valid. 

Moving average smoothing is often thought of as 
a crude low pass digital filter [Ref. 9]. The smoothing 
operation filters out the high frequency’ turbulent 
fluctuations leaving a low frequency "local" mean. As the 
length of the moving average window is increased (increasing 
L), lower and lower frequencies are filtered out. So, the 
resulting mean velocity based on the ensemble average of the 
smoothed velocity is partly dependent on the frequency 
response characteristics of the smoothing filter. The problem 
becomes one of choosing an appropriate smoothing window size 
to obtain a specific low pass filter cutoff frequency. It 
will be left up to the aerodynamicist to choose the 
appropriate cutoff frequency based on knowledge of the nature 


of unsteady turbulent boundary layers. However, the 
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frequencies of interest may be bounded. Since the propeller 
is generating turbulent pulses at a frequency of approximately 
40 Hz, the cutoff frequency should be greater than 40 Hz. 
Landrum and Macha [Ref. 14] show that boundary layers in 
transition from laminar to turbulent flow contain spectral 
peaks at approximately 800 Hz. So, the cutoff frequency of 
a low pass filter should be less than 800 Hz to filter out the 
fluctuating component. Once a cutoff frequency is selected, 
all subsequent analysis should be accomplished using only the 
selected value. 

Table IV-1 presents the relationship between the 
size of the moving average window and the cutoff frequency of 
the moving average digital filter for a sample rate of 15,000 
and 30,000 Hz. The cutoff frequency is defined as the 
frequency at which the magnitude of the filtered signal has 
been attenuated by -3 dB (a factor of 0.707) relative to the 
unfiltered signal. Included in Appendix A is a derivation and 
discussion of the frequency response characteristics of the 
moving average digital filter. As detailed in Appendix A, the 
moving average smoothing algorithm possesses poor digital 
filter characteristics including a severe Gibb's phenomenon 
(positive and negative ripples in the frequency response 


function, see Appendix A). 
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Table IV-1. APPROXIMATE CUTOFF FREQUENCIES FOR MOVING 
AVERAGE WINDOWS 


Window Cutoff Freq. (Hz) Cutoff Freq. (HZ) 
Size (Sample rate=15000 Hz) (Sample rate=30000 Hz) 





The moving average algorithm for determining the 


non-stationary mean velocity 1s accomplished through the 


following steps: 


Identify the input and output data file names and enter 
the desired moving average window size (an odd number of 
points). 


Separate the input time history record into an ensemble 
of individual pulses. As previously described, each 
column of array ENSMBL represents a single propeller 
pulse. 


Apply moving average smoothing to each column (propeller 
pulse) of array ensemble. Write the smoothed values of 
velocity back into array ENSMBL for storage. 


Perform ensemble averaging on the smoothed velocity to 
obtain the final mean velocity. 


Write the output file. 


.e) Il 


b. Implementation and Examples 


The moving average algorithm was implemented as 
program MOVAVE. The source code and user's guide for program 
MOVAVE is presented in Appendix B. Using test case 1 as the 
input file, program MOVAVE produced the results presented in 
Figure IV-6. A 21 point moving average window was selected 
yielding a cutoff frequency of 320 Hz. Comparing the moving 
average results of Figure IV-6 to the pure ensemble average 
shown in Figure IV-5, observe that, as expected, the small 
high frequency fluctuations remaining in the ensemble average 
have been removed. Figure IV-7 presents the results of 
program MOVAVE with a 7 point smoothing window selected (960 
Hz. cutoff frequency). Observe that some high frequency 
Signal is still present in the mean velocity suggesting that 
a Significant amount of turbulent energy exists below 960 Hz. 
(Later spectral analysis will show this to be true.) Figure 
IV-8 shows the results of program MOVAVE with a 41 point 
moving average window (corresponding to a cutoff frequency of 
1G Siete ee Note that this example is not significantly 
different from the 21 point window results shown in Figure IV- 
6. The similarity of Figures IV-6 and IV-8 suggests that the 
21 point window has succeeded in removing most of the high 
frequency fluctuating velocity component. 

Combining ensemble averaging with a "low 


frequency mean" empirical model (moving average smoothing) 
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Figure IV-6. RESULT OF PROGRAM MOVAVE ON TEST CASE 1. 
(SMOOTHED USING A 21 POINT MOVING AVERAGE, THEN ENSEMBLE 


AVERAGED. ) 


800 


T CVT Tas 


600 







—- 960 Hz 


Cutoff Freq. 


(Pe ee ee eee ee ee, 
400 
Sample Number 


200 


© © Oo 
© © 
wv 


WN 
(syunod. [eyIdIp) AzIOo[aA, ueapy 


600 





Figure IV-7. RESULT OF PROGRAM MOVAVE ON TEST CASE 1. 
(SMOOTHED USING A 7 POINT MOVING AVERAGE, THEN ENSEMBLE 
AVERAGED. ) 
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Figure IV-8. RESULT OF PROGRAM MOVAVE ON TEST CASE 1. 
(SMOOTHED USING A 41 POINT MOVING AVERAGE, THEN ENSEMBLE 


AVERAGED. ) 
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provides at least a qualitative improvement in the estimate 
of the non-stationary mean velocity without increasing the 
volume of data. Since the mean of the non-stationary velocity 
is somewhat arbitrarily defined as the “low frequency part", 
it becomes difficult to quantify any improvement in the 
accuracy of the estimate of the mean. The resulting mean 
velocity is strongly dependent upon the size of the moving 
average smoothing window (i.e. the cutoff frequency). The 
cutoff frequency must be carefully selected to be consistent 
with the low frequency mean empirical model and the physical 
nature of the unsteady turbulent boundary layer. 
4. Frequency Domain Smoothing Algorithm 


a. Description of the Algorithm 


Frequency domain smoothing involves the more 
classical approach to digital filtering. The time history 
record is transformed into the frequency domain using the 
discrete Fourier transform. The frequency response function 
of the desired low pass filter is multiplied by the 
transformed data to remove the high frequency content. Then 
the data are transformed back into the time domain as the 
smoothed (low pass filtered) velocity. Recall that 
multiplication in the frequency domain is equivalent to 
convolution in the time domain. 

The ideal low pass filter (in the frequency 


domain) looks like Figure IV-9, where the step occurs at the 
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desired cutoff frequency (f,). As discussed in section 
III.B.6., practical finite length time series cause errors in 
the Fourier series representations of functions due to 
truncation of the Fourier coefficients. Therefore the "ideal" 
filter cannot be achieved practically. Figure IV-10 
illustrates the result of the truncated Fourier series 
representation of the ideal filter. The ripples in the 
frequency response function are known as the Gibb's phenomenon 
(Ref. 9]. The Gibb's phenomenon 1S an_ undesirable 
characteristic for a filter. Smooth, continuous filter 
frequency response functions are selected to minimize the 
ripple effect (i.e. to reduce the magnitude of the high 


frequency Fourier coefficients per section III.B.6.). 
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Figure IV-9. THE IDEAL LOW PASS FILTER. 

The specific filter frequency response 
function used in the smoothing algorithm presented here, is 
a parabolic shaped function. The filter function is one at 
zero frequency and falls smoothly to zero at a frequency which 


depends on the user specified smoothing window size (the 
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Figure IV-10. TRUNCATED FOURIER SERIES REPRESENTATION OF AN 
IDEAL FILTER. 


cutoff frequency). An analysis of the frequency response 
characteristics of the parabolic filter is presented in 
Appendix A. Other frequency domain filters with different 
frequency response characteristics may be implemented, such 
as the von Hann window described in section III.B.6. Table 
IV-2 presents the cutoff frequency, (f.), as a function of the 
smoothing window size, (L), the sample rate, (f,), and the 
size (number of points) of the DFT, (M). 

The frequency domain smoothing algorithm for 
determining the non-stationary mean velocity is very similar 
to the moving average algorithm. Frequency domain low pass 
filter smoothing is accomplished prior to ensemble averaging 
instead of a moving average. The low pass filter algorithm 
requires a program to perform the DFT and IDFT. The following 
steps describe the algorithm: 


- Identify the input and output data file names and enter 
the desired frequency domain smoothing window size. 
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Table IV-2. APPROXIMATE CUTOFF FREQUENCIES FOR FREQUENCY 
DOMAIN LOW PASS FILTER. 


tee (Hi2") f° (Hz) 
f==15000, M=1024) (f5=30000, M=2048) 


2 
3 
4 
5 
6 
7 
8 
=) 





* Separate the input time history record into an ensemble 


of individual pulses. As previously described, each 
column of array ENSMBL represents a single propeller 
pulse. 


- Apply the low pass filter smoothing to each column 
(propeller pulse) of array ensemble. This step requires 
an FFT routine to accomplish the DFT and IDFT. Write the 
smoothed values of velocity back into array ENSMBL for 
storage. The smoothing subroutine and the FFT subroutine 
were extracted from Press, et al. [Ref. 12]. 


- Perform ensemble averaging on the smoothed velocity to 
obtain the final mean velocity. 


¢ Write the output file. 
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b. Implementation and Examples 


The frequency domain smoothing algorithm was 
implemented as program SMUMEAN. The source code and user's 
guide for program SMUMEAN is included in Appendix B. The low 
pass filter and FFT subroutines used in this implementation 
are from Press, et al. [Ref. 12]. Figure IV-11 presents the 
results of applying program SMUMEAN on test case 1 for a 25 
point smoothing window (a cutoff frequency of 325 Hz). 
Figures IV-12 and IV-13 present the results of program SMUMEAN 
for an 8 point (f, = 1015 Hz) and a 51 point (f, = 160 Haz) 
smoothing window respectively. These cutoff frequencies were 
chosen as close as possible to those presented in Figures IV- 
6 through iV-s. 

Comparisons of the frequency domain smoothing 
results with the moving average smoothing results indicate 
negligible differences in the corresponding mean velocity 
between the two algorithms for similar cutoff frequencies. 
These results are consistent with the fact that the frequency 
response function of the frequency domain smoothing algorithm 
is within 1% of the frequency response function of the moving 
average algorithm at frequencies between zero and the cutoff 
frequency (see Appendix A). The frequency response functions 
of the two methods begin to diverge significantly above the 
cutoff frequency. But at this point, the magnitude of the 


high frequency fluctuating velocity component has been 
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Figure IV-11. RESULT OF PROGRAM SMUMEAN ON TEST CASE 1. 
(SMOOTHED USING A 25 POINT LOW PASS FILTER, THEN ENSEMBLE 
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Figure IV-12. RESULT OF PROGRAM SMUMEAN ON TEST CASE 1. 
(SMOOTHED USING AN 8 POINT LOW PASS FILTER, THEN ENSEMBLE 


AVERAGED. ) 
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Figure IV-13. RESULT OF PROGRAM SMUMEAN ON TEST CASE 1. 
(SMOOTHED USING A 51 POINT LOW PASS FILTER, THEN ENSEMBLE 
AVERAGED. ) 
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significantly attenuated such that there is little effect on 
the resulting mean velocity. Choosing a low pass filter for 
the frequency domain smoothing algorithm with a frequency 
response function that better approximates the "ideal" low 
pass filter would cause a greater (though still minor) 
difference between the two algorithms. 

Combining ensemble averaging with frequency domain 
smoothing (low pass filtering) yields similar results compared 
to the moving average smoothing algorithm. A similar 
improvement is obtained in the estimate of the mean of a non- 
stationary velocity measurement as was obtained with program 
MOVAVE, when compared to ensemble averaging alone. Again, the 
cutoff frequency must be chosen with care because the 
resulting mean velocity is strongly dependent upon the cutoff 
frequency selected. The frequency domain smoothing algorithm 
has the advantage of increased flexibility in the choice of 
the low pass filter characteristics. Furthermore, frequency 
domain smoothing has a larger choice of cutoff frequencies 
(compare Table IV-1 with Table IV-2), and the beginning and 
end points are not lost (set to zero) as occurs using the 
moving average algorithm. 

5. Ensemble Averaging Before Smoothing 
Both of the algorithms described in sections IV.B.3. 
and IV.B.4 above employ smoothing (low pass’ filtering) 


techniques prior to ensemble averaging. Both moving average 
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smoothing and frequency domain smoothing are linear 
operations, as is ensemble averaging. Therefore, identical 
mean velocity results are obtained when the ensemble average 
operation is performed prior to the smoothing operation. 
However, a substantial increase in computational efficiency 
is achieved because the smoothing operation is only performed 
on a single ensemble averaged velocity record instead of on 
each of the individual pulses which comprise the ensemble. 

Program MEANSMU was implemented to demonstrate the 
increased computational speed achieved by ensemble averaging 
prior to smoothing. The algorithm for program MEANSMU is the 
same aS program SMUMEAN presented in section IV.B.4. above, 
except that steps three and four are reversed. Appendix B 
presents the source code and user's guide for program MEANSMU. 
The results presented in Figure IV-11 required approximately 
143 seconds of computer time for program SMUMEAN. Program 
MEANSMU produces identical results in only 34 seconds of 
computer time. (Note that the computer was a ten megahertz 
IBM PC/AT clone equipped with an 80287 math co-processor). 
No results are presented for program MEANSMU since they are 
identical to those shown in Figures IV-11 through IV-13. 

6. Summary--Mean Velocity Determination 

Program ENSEMBL provides a method for estimating the 

non-stationary mean velocity of an ensemble of propeller blade 


wake passages by ensemble averaging. Program MOVAVE and 


Ves, 


SMUMEAN provide alternate methods for obtaining improved 
estimates of the non-stationary mean velocity when the number 
of pulses (i.e. the number of sample records) in the ensemble 
is small. Program MOVAVE and SMUMEAN yield similar results. 
Program SMUMEAN maintains the following advantages over 
program MOVAVE: 


+ Increased flexibility in the number of available cutoff 
frequencies. 


* Increased flexibility in the choice of the frequency 
response characteristics of the low pass filter used for 
smoothing. 


- No loss of data at the beginning and end of the mean 
velocity data record. 


Program MEANSMU possesses these same advantages plus the added 
advantage of increased computational efficiency. Therefore, 
program MEANSMU is the recommended method for determining the 
non-stationary mean velocity. It is further recommended that 
alternate frequency domain low pass filters be implemented and 
tested in the "smoothing" algorithm of program MEANSMU to 
determine if there is any significant difference or 


improvement in the non-stationary mean velocity estimate. 


C. TURBULENCE INTENSITY DETERMINATION 


1. Overview 
In Chapter II, turbulence intensity was defined as 
the root mean square of the fluctuating velocity component. 


For steady turbulence, the fluctuating velocity component is 
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simply the difference between the instantaneous velocity and 
the mean velocity. Thus, for steady turbulence, turbulence 
intensity is very similar to the statistical variance. But 
for unsteady (non-stationary) turbulence, the fluctuating 
velocity component is defined as the difference between the 
instantaneous velocity and the local smoothed velocity (the 
local mean), not the non-stationary ensemble averaged mean 
velocity. The local smoothed velocity is obtained using 
either the moving average smoothing or the frequency domain 
low pass filter smoothing. And the "mean" operation in the 
root mean square of the fluctuating component is the ensemble 
average. This definition prevents variations in the local 
mean from pulse to pulse to be incorrectly included in the 
turbulence intensity measure. In other words, for unsteady 
turbulence, turbulence intensity measures only the variation 
due to the fluctuating velocity component relative to the 
smoothed velocity (local mean). But statistical variance 
measures the total variation of the instantaneous velocity 
due to variations in both the fluctuating velocity component 
(turbulence) and variations in the local smoothed velocity. 
Figure IV-14 presents an illustration. 

Using this "local mean" definition of turbulence 
intensity, two algorithms have been developed using the two 


smoothing algorithms which were presented in section IV.A. 
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Figure IV-14. ILLUSTRATION OF TURBULENCE INTENSITY (BASED ON 
"LOCAL MEAN'') VERSUS STATISTICAL VARIANCE (BASED ON ENSEMBLE 
AVERAGE OF LOCAL MEAN) 


above. The turbulence intensity algorithms are presented in 
the following paragraphs. 
2. Turbulence Intensity from Moving Average Smoothing 


a. Description of the Algorithm 


Moving average smoothing was described in section 
IV.A.3. above. The local mean (smoothed velocity) at a given 


gleheayore ale Veabw le UlAes) is used to determine the local 


LOCAL / 
fluctuating velocity component. Expanding on equation II-3, 


turbulence intensity is computed as: 
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ix . 2 

Re u(t) ate ise 

TURBULENCE INTENSITY = - (IV-1) 
edge 


where k is the ensemble index, N is the total number of 
propeller pulses, and U,,, 1s a reference velocity (usually 
the mean velocity at the edge of the boundary layer) used to 
normalize the turbulence intensity. U ace is constant for a 
particular boundary layer profile. 

The assumptions which apply to the moving average 
mean velocity algorithm also apply to the moving average 
turbulence intensity algorithm. The instantaneous velocity 
data is assumed stationary over the short smoothing window 
such that time averaging over the short segment may be used 
to estimate the local mean. Also the data are assumed to 
behave according to the "low frequency mean" empirical model. 
Therefore, as with the moving average mean velocity algorithn, 
the turbulence intensity results are a strong function of 
cutoff frequency (smoothing window size). Table IV-1 presents 
the relationship between moving average window size and cutoff 
frequency. As before, the choice of an appropriate cutoff 
frequency will be left up to the aerodynamicist. The 
frequency response characteristics of moving average smoothing 


is presented in Appendix A. For consistent results, the same 
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cutoff frequency should be used in the turbulence intensity 
analysis as was used in the mean velocity analysis. 
The moving average turbulence intensity algorithm 
is accomplished through the following steps: 
- Identify the input and output data file names and enter 
the desired moving average window size (an odd number of 
points). 


¢ Separate the input time history record into an ensemble 


of individual pulses. As previously described, each 
column of array ENSMBL represents a single propeller 
pulse. 


¢ Apply moving average smoothing to each column (propeller 
pulse) of array ensemble. Subtract the instantaneous 
velocity from the smoothed velocity to get the 
fluctuating velocity component. Square the fluctuating 
velocity component and write the result back into array 
ENSMBL for storage. 

- Perform ensemble averaging on the squared fluctuating 
velocity to obtain the mean square fluctuating velocity. 
Then compute the square root and normalize (divide by 
Usdge) 0 obtain the turbulence intensity. 

- Write the output file. 


b. Implementation and Examples 


The moving average turbulence intensity algorithm 
was implemented as program TURINTMA. The source code and 
user's guide for program TURINTMA is presented in Appendix B. 
Figures IV-15, IV-16, and IV-17 present the results of program 
TURINTMA (using test case 1) for cutoff frequencies of 320 Hz, 
960 Hz and 165 Hz respectively. These cutoff frequencies 
correspond to moving average window sizes of 21 points, 7 


points and 41 points respectively. The turbulence intensity 
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results presented in Figures IV-15, IV-16 and IV-17 correspond 
to the mean velocity results presented in Figures IV-6, IV-7 
and IV-8. These data show the time dependent, periodic nature 
of the unsteady turbulent boundary layer as the fluid flow 
oscillates between the turbulent and laminar states. 
In comparing the turbulence intensities of Figures 
IV-15, IV-16 and IV-17, observe that as cutoff frequency is 
decreased, turbulence intensity level increases. This 
illustrates the strong relationship between cutoff frequency 
and turbulence intensity. The cutoff frequency must be 
carefully selected to be consistent with the low frequency 
mean empirical model and the physical nature of the unsteady 
turbulent boundary layer. Note that the moving average 
smoothing has caused the loss of a small segment of data at 
the beginning and end of the data record. 
3. Turbulence Intensity from Frequency Domain Smoothing 


a. Description of Algorithm 


Frequency domain low pass filter smoothing was 
described in section IV.A.4. above. The turbulence intensity 
algorithm using frequency domain smoothing is essentially the 
same as that using moving average smoothing. The key 
difference is the use of low pass frequency domain filtering 
to obtain the smoothed (local mean) velocity, u(t;),o4,- The 
turbulence intensity is still computed using equation IV-1. 


Table IV-2 presents the relation between the smoothing window 
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Figure IV-15. RESULT OF PROGRAM TURINTMA ON TEST CASE 1. 
(TURBULENCE INTENSITY COMPUTED USING 21 POINT MOVING AVERAGE 


SMOOTHING. ) 
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Figure IV-16. RESULT OF PROGRAM TURINTMA ON TEST CASE 1. 
(TURBULENCE INTENSITY COMPUTED USING 7 POINT MOVING AVERAGE 
SMOOTHING. ) 
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Figure IV-17. RESULT OF PROGRAM TURINTMA ON TEST CASE 1. 
(TURBULENCE INTENSITY COMPUTED USING 41 POINT MOVING AVERAGE 


SMOOTHING. ) 
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size and the cutoff frequency. The frequency response 
characteristics of the low pass filter used in this algorithm 
are presented in Appendix A. 

The algorithm outlined below describes the steps 
required to determine the turbulence intensity using frequency 
domain low pass filter smoothing: 

- Identify the input and output data file names and enter 
the desired low pass filter window size corresponding to 
the desired cutoff frequency of the filter. 

¢ Apply low pass filter smoothing to each column (propeller 
pulse) of array ensemble. This step requires an FFT 
routine to accomplish the DFT and IDFT. Subtract the 
instantaneous velocity from the smoothed velocity to get 
the fluctuating velocity component. Square the 
fluctuating velocity component and write the result back 
into array ENSMBL for storage. 


* Perform ensemble averaging on the squared fluctuating 
velocity to obtain the mean square fluctuating velocity. 


* Then compute the square root and normalize (divide by 
Usgge) tO obtain the turbulence intensity. 


¢ Write the output file. 


b. Implementation and Examples 


The algorithm for computing turbulence intensity 
using frequency domain smoothing was implemented as program 
TURBIN. The source code and user's guide for program TURBIN 
is included in Appendix B. The low pass filter and FFT 
subroutines used in this implementation are from Press, et al. 
Peet. 12}. Figure IV-18 presents the results of applying 
program TURBIN to test case 1 for a 25 point smoothing window 


(a cutoff frequency of 325 Hz). Figures IV-19 and IV-20 
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present the results of program TURBIN for an 8 point ss@je 
1015 Hz) and a 51 point (f, = 160 Hz) smoothing window 
respectively. These cutoff frequencies were chosen as close 
aS possible to those presented in Figures IV-15 through IV- 
17. The turbulence intensity results shown in Figures IV-18 
through IV-20 correspond to the mean velocity results shown 
in Figures IV-11 through IV-13 respectively. Again, the 
cutoff frequency must be chosen with care because the 
resulting turbulence intensities are strongly dependent upon 
the cutoff frequency selected. 

Since the frequency response characteristics of 
moving average smoothing are similar to those of frequency 
domain smoothing (for the particular filter implemented here), 
the turbulence intensity results are also very similar. 
Respective comparison of Figures IV-18, IV-19 and IV-20 to 
Figures IV-15, IV-16 and IV-17 illustrate the similarity. 
Analogous to the mean velocity algorithms, turbulence 
intensity from frequency domain smoothing has the advantage 
of increased flexibility in the choice of the low pass filter 
characteristics compared to moving average smoothing. Also, 
frequency domain smoothing has a larger choice of cutoff 
frequencies, and the beginning and end points are not lost 


(set to zero) as occurs using the moving average algorithm. 
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Figure IV-18. RESULT OF PROGRAM TURBIN ON TEST CASE 1. 
(TURBULENCE INTENSITY COMPUTED USING 25 POINT FREQUENCY DOMAIN 
LOW PASS FILTERING.) 
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Figure IV-19. RESULT OF PROGRAM TURBIN ON TEST CASE 1. 
(TURBULENCE INTENSITY COMPUTED USING 8 POINT FREQUENCY DOMAIN 
LOW PASS FILTERING. ) 
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Figure IV-20. RESULT OF PROGRAM TURBIN ON TEST CASE 1. 
(TURBULENCE INTENSITY COMPUTED USING 51 POINT FREQUENCY DOMAIN 
LOW PASS FILTERING.) 
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4. Summary--Turbulence Intensity Determination 

Programs TURINTMA and TURBIN provide two alternate 
methods for computing the non-stationary turbulence intensity 
of unsteady turbulent boundary layer velocity data. AS 
implemented, the two programs produce nearly identical results 
for the same cutoff frequency. The advantages of the low pass 
filter smoothing technique are described in the paragraph 
immediately above. These advantages lead to the 
recommendation that program TURBIN be selected as the primary 
method for determination of the non-stationary turbulence 
intensity. It is also recommended that alternate frequency 
domain low pass filters be implemented and tested in the 
"smoothing" algorithm of program TURBIN to determine if there 
is any Significant difference or improvement in the non- 
stationary turbulence intensity results. In particular, 
filters which better approximate the "ideal" low pass filter 


should be examined. 


D. SPECTRAL ANALYSIS 


1. Overview 
The primary purpose of this section 1s to describe 
the technique developed to estimate the time varying frequency 
spectra of non-stationary data. Spectral analysis algorithms 
for stationary data are developed initially to gain an 


understanding of the techniques involved. Two definitions of 
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the spectral density function were introduced in Chapter III. 
The first defined the power spectral density (PSD) function 
as the Fourier transform of the autocorrelation function. The 
second, more direct definition obtains the PSD by computing 
the squared magnitude of the Fourier transform of the time 
history record. The second definition is the one developed 
in this section. Averaging techniques are employed to obtain 
a stable, minimum variance spectral estimator. 
2. Spectral Analysis of Stationary Data 


a. Description of the Algorithm 


The spectral estimation scheme developed here is 
based on the classical "overlapping segments" method developed 
by Welch and described in Marple [Ref. 10}. The Welch method 
was selected because it provides an efficient FFT based 
algorithm. Also, the Welch method does not depend on the 
assumption of a specific model for the input data (the Welch 
method is "non-parametric" [Ref. 10]. 

The idea is to divide the time history record into 
a number of shorter segments, as illustrated in Figure IV-21. 
Each segment is multiplied by a window function of the same 
length as the segment. The purpose of the tapering window 
function is to minimize the frequency "leakage" effect of the 
Gibb's phenomenon caused by truncation of the Fourier series 
(as outlined in section III.B.6.). The "raw" PSD estimate of 


each windowed segment is then computed per equation III-35. 


ou 


The power in each segment is corrected for the window function 
bias by dividing each segment PSD by the window function 
energy (the sum of the square of each window function term). 
The PSD estimates for all the segments are then averaged to 
produce the final Welch PSD estimate. The spectrum is 
computed from a frequency of zero to one half the sample rate 


(i.e., the Nyquist frequency, see section III.B.2). 
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Figure IV-21. SEQUENCE OF OVERLAPPED SEGMENTS OF LENGTH WITH 
50% OVERLAPPING. 





This "segment and average" method of PSD 
estimation is necessary because it reduces the variance of the 
PSD estimate. A raw PSD estimate (equation III-35) has been 
shown to have a variance on the same order as the estimate 
itself, making the estimate unreliable at best [Ref. 10]. 
Marple shows that the variance of the Welch PSD estimator is 
roughly inversely proportional to the number of segments. So 
ideally, a large number of segments is desirable. However, 
the tradeoff is that using a large number of segments reduces 


the frequency resolution (the number of discrete frequency 


iz 


bins). The number of positive frequency bins is equal to half 
the number of data points in the segment. Thus the segments 
cannot be made too small. For 50% overlapping segments, the 
relation between the number of positive frequency bins, M, the 
number of segments, 2K, and the total number of data points 
in the time series, N, is given by the equation: 

N = (2K +1)M (IV-2) 
Since an FFT algorithm is used to compute the DFTs, M must be 
an integer power of two and N ideally should be an integer 
power of two. (Note that due to overlapping, some values of 
K and M yield a value of N which is slightly less than an 
integer power of 2. The required number of input data points, 
N, for the spectral analysis algorithm will never be greater 
than that specified in Table IV-3.) Table IV-3 presents some 
typical values of the parameters N, M and K. The goal is to 
choose an adequate number of frequency bins, M with the 
largest number of segments, K possible. At least 14 segments 
(K=7) is desirable. The easiest solution to the tradeoff is 
to collect longer data records. Then both a sufficiently 
large K and M may be specified. 

Many different window functions have _ been 
developed to reduce the spectral leakage phenomenon described 
in Chapter III. Each window function possesses different 
frequency response characteristics. Window functions are 


selected based on a tradeoff between frequency resolution loss 
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Table IV-3. TYPICAL PARAMETER VALUES FOR THE WELCH PSD 
ESTIMATOR 





(main lobe width) and the amount of leakage suppression 
(highest side lobe level). An outstanding comparison and 
discussion of the characteristics and relative merits of most 
of the popular window functions may be found in the paper by 
Harris (Ref. LBS al lpe Seven of these window functions were 
selected to be incorporated into the spectral estimation 
algorithm and may be selected by the user. Table IV-4 lists 
the specific window functions with some of their frequency 
response parameters. The Hamming window or the three term 
Blackman-Harris window provide a good compromise. The 
formulas for the window functions may be found in Marple [Ref. 
10] or Harris [Ref. 13] and are presented in Appendix B in the 


FORTRAN function subroutine called FUNCTION WINDOW. 
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Table IV-4. WINDOW FUNCTIONS 


Highest Side 6 dB 
Lobe Level Bandwidth 
Window (dB) (bins) 


Parzen (triangular) 


Rectangular (no window) 
Welch (eh se Jere) Lae) 

van Hann (cosine“) 
Hamming 

3 term Blackman-Harris 
4 term Blackman-Harris 





The Welch spectral estimation algorithm is 


accomplished through the steps listed below: 


e 


Enter the input data file name, the output data file name 
and the sample rate of the data. Select the window 
function option and enter the parameters for the number 
of segments (2K), and the number of frequency bins (M) 
based on the integer power of two number of data points 
(N). The input data file must contain at least n data 
points where N is computed from equation IV-2. 


Allow selection of the option to remove the stationary 
mean value from the data (i.e., the constant mean 
velocity of a steady turbulence measurement). Data with 
a large mean value will yield a spectral estimate with 
a very large power magnitude at zero frequency. Due to 
leakage, which windowing reduces but does not eliminate, 
the large zero frequency power may mask the true spectral 
characteristics of nearby (low) frequencies. If this 
option is selected, compute the mean velocity and 
subtract the mean from each of the instantaneous velocity 
values. Write the resulting fluctuating velocity 
component to a scratch file to be used for spectral 
estimation. 


Compute the spectral estimate using the Welch 
"overlapping segments" method. This is accomplished 
using a modified version of a FORTRAN’ subroutine 
(SUBROUTINE SPECTRM) extracted from Press et al. [Ref. 
12]. The FFT subroutine was also obtained from Reference 
a2. 
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* Compute the relative power in dB, of the resulting PSD 
estimate using the peak value as the reference (zero dB) 
point. Also compute the percent running total power at 
each increasing frequency relative to the total power 
summed over all frequencies. 

- Write the output file 


b. Implementation and Examples 


The spectral estimation algorithm for stationary 
data described above was implemented as program PSD. The 
source code and user's guide for program PSD is presented in 
Appendix B. Figure IV-22 presents the results of program PSD 
on test case 1. Though test case 1 is not strictly 
stationary, program PSD will provide a "time averaged" 
spectral estimate over the length of the time history record. 
The Hamming spectral window was selected for this example. 
Program PSD was executed with 256 frequency bins selected 
(M=256) and 62 segments selected (K=31). These parameter 
values correspond to 16,128 data points or slightly more than 
half of the total number of points in the time history record. 
To obtain a PSD estimate of the entire time history record, 
the input file may be "zero padded" to a total length equal 
to an integer power of two data points (32,768 in this case). 
That is, add zeros on the end of the input file to obtain the 
proper "power of two" length. 

The option to remove the mean was selected for 
this example (mean value = 317 digital velocity counts) 


because the resulting zero frequency power spectrum value 
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Figure IV-22. SPECTRAL ESTIMATE OF TEST CASE 1 TREATED AS A 
STATIONARY SIGNAL. (HAMMING WINDOW, STATIONARY MEAN REMOVED, 
FREQUENCY BINS--256, OVERLAPPED SEGMENTS--31, SAMPLE RATE-- 
15,000 HZ.) 
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would otherwise be extremely large (on the order of 10) 
compared to the remainder of the frequency spectrum. Even 
with the mean removed, the power spectrum contains a high peak 
centered at the propeller pulse rotational frequency of 40 Hz. 
The interesting part of the frequency spectrum is well above 
the propeller pulse frequencies where the turbulent 
fluctuations lie. Therefore, the spectrum data below 300 Hz. 
was not shown and the vertical scale was enlarged to better 
illustrate the high frequency portion of the spectrum. An 
alternate approach would be to use the low pass filtering 
algorithm (frequency domain smoothing) to compute and then 
remove the low frequency velocity component. This would be 
equivalent to high pass filtering. Finally, program PSD is 
used without mean removal (since the mean has already been 
removed) to estimate the spectrum. 

Test case 1 was also analyzed using the three term 
Blackman-Harris window which has lower peak sidelobes than the 
Hamming window. The resulting spectral estimate plot is not 
shown because it is virtually identical to Figure IV-22. The 
three term Blackman-Harris window may be better suited to 
observing small (low power) spectral peaks which are widely 
separated from larger peaks. However, Figure IV-22 
1llustrates the apparent broad band nature of the boundary 


layer turbulence. Therefore, any of the popular window 
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functions (options four through seven in Table IV-4) will 
produce much the same result. 


3. Spectral Analysis of Stationary Data--Maximum Entropy 
Method 


An alternate, parametric spectral estimation algorithm 
was briefly examined. Called the "all poles" or "maximum 
entropy" method (MEM), this algorithm is based on an 
autoregressive (AR) parametric model of the time history data 
(Ref. 10]. The resulting spectral estimates produce spectral 
peaks which are no longer linearly proportional to the power 
contained in the underlying sinusoid. Therefore, MEM spectral 
estimates cannot be averaged as was done using classical 
methods. The shear volume of data involved with unsteady 
turbulent boundary layer analysis require the use of 
averaging. Therefore, MEM spectral estimation cannot be used 
for the non-stationary application which will be developed in 
the next section. 

Program PSDMEM is presented in Appendix B for the sake 
of completeness. The program is included "as is" with no 
detailed explanations or documentation. This program may be 
used to obtain spectral estimates of stationary data. The 
maximum entropy method is especially suited to extracting 
narrow spectral peaks which are close together. (Which is not 
important to turbulence study.) Also MEM spectral estimation 
does not require an integer power of two data points and can 


provide reasonable estimates with a relatively small number 
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of data points. The program was extracted from Press et al. 
(Ref. 12]. Additional theory on the maximum entropy method 
may be found in Marple [{Ref. 10}. 

4. Spectral Analysis of Non-Stationary Data 


a. Description of the Algorithm 


With non-stationary data, the frequency content 
of the data varies with time. For the particular application 
of unsteady periodic boundary layer velocity data, ensemble 
averaging techniques may be employed on an ensemble of 
propeller "pulses". The algorithm is somewhat analogous to 
the moving average scheme. After the data are separated into 
an ensemble of individual pulses, a raw PSD estimate is 
computed for a short time segment from the first pulse per 
equation III-35. This represents the spectrum of the first 
pulse at the center time point of the segment. Prior to 
computing the raw spectral estimate, the segment is windowed 
using an appropriate tapering window just as in the algorithm 
for stationary data above. Similarly, a raw PSD estimate is 
computed for the corresponding time segment in every other 
pulse of the ensemble, windowing the data as before. The 
collection or raw PSD estimates from the ensemble of short 
time segments is then averaged and normalized to produce the 
final "instantaneous" spectral. estimate for the ensemble at 


the time represented by the center of the segment. As before, 
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the purpose of the normalization is to correct the power for 
the effect of the tapering window. 

The above process can then be repeated for other 
ensembles of short time segments producing ensemble averaged 
spectral estimates centered at different times across the 
ensemble. In this manner, the time varying spectral 
characteristics may be constructed. The assumption associated 
with this algorithm (as with the moving average algorithm) is 
that the data are stationary over the short time segment. 
This assumption is reasonable since the mean varies slowly 
compared to the high frequency fluctuating velocity component. 
However, the time segment must be long enough to achieve the 
desired resolution in the frequency spectrum. If the 
frequency content of the data does change in time over the 
segment, then the spectral estimate is still useful. The 
spectral estimate then represents the average frequency 
content of the data over the time segment rather than the 
instantaneous spectrum at the center of the time segment. 

Evaluating the tradeoff between short segment 
length and desired frequency resolution led to the choice of 
256 point segments. This provides 128 positive frequencies 
in the frequency spectrum which is adequate. Also, 256 points 
represents roughly 1/3 to 1/6 the length of the entire pulse 
for sample rates of 15,000 and 30,000 Hz. respectively. A 


segment which is 1/3 the length of the entire pulse is 
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considered marginal for the "stationary short segment" 
assumption. Increasing the sample rate to 20,000 Hz. or more 
is an easy solution to the problem. Another equivalent 
alternative is to decrease the frequency resolution to 64 
positive frequencies which corresponds to a segment width of 
128 points. Yet another option is to maintain a 256 point 
wide window, but (similar to the moving average window) 
increase the resolution in time by overlapping the ensemble 
of segments such that the final spectral estimates are 
centered on time points which are less than 256 points apart 
(say, 50 or 100 points for example). These options are not 
mutually exclusive and may be employed together. 

The variance of the final spectral estimate of a 
segment of the ensemble is inversely proportional to the 
number of pulses which comprise the ensemble. This is 
analogous to the stationary case except the segments are in 
the ensemble sense and not overlapping segments in time. 
Therefore 30 to 50 pulses should be adequate to obtain a 
reasonable stable spectral estimate. This would roughly 
correspond to a K parameter in the range of 15 to 25. 

This non-stationary spectral estimation algorithm 
provides the same seven tapering window options which were 
described above for the stationary case. The window functions 
and some of their frequency response characteristics are 


presented inn Table IV-4. As before, the Hamming window or 
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the three term Blackman-Harris window are recommended for the 


broad band turbulence data. 


The non-stationary spectral estimation algorithm 


is summarized using the following steps listed below: 


Enter the input data file name, the output data file name 
and the sample rate of the data. 


Separate the input time history record into an ensemble 
of individual pulses. As previously described, each 
column of array ENSMBL represents a single propeller 
pulse. 


Allow selection of the option to remove the mean value 


from the data (i.e., remove the mean velocity of each 
individual segment). Data with a large mean value will 
yield a spectral estimate with a very large zero 
frequency value. This may mask the true spectral 


estimate of the nearby (low) frequencies. 
Allow selection of the tapering window option. 


Allow selection of one of three options for segmenting 
the ensemble. The first option segments the ensemble 
based on an input value for the number of data points 
between the beginning of each segment (i.e., the amount 
of time between each PSD). The second option segments 
the data based on an input value for the number of 
segments across the ensemble of pulses (1.e., the number 
of PSDs across the ensemble). The third option creates 
a Single segment ensemble based on an input time (sample) 
value of the center of the segment (1i.e., choose the 
location for a single PSD). 


Remove the mean of each segment if selected above. Then, 
compute the raw spectral estimate of the corresponding 
segment of each pulse, windowing the data prior to 
calling the FFT subroutine. The FFT subroutine was 
extracted from Press, et al. [Ref. 12]. Now average the 
raw PSDs from each pulse and normalize to correct for the 
window. 


Compute the relative power of the final segment PSD in 


dB using the peak value as the reference (zero dB) point. 
Also compute the percent running total power at each 


a O3 


increasing frequency relative to the total power summed 
over all frequencies. 
¢ Write the output file. 


b. Implementation and Examples 


The algorithm for non-stationary spectral 
estimation was implemented as program ENSPSDAV. The source 
code and user's guide is presented in Appendix B. Figure IV- 
23 presents the results of applying program ENSPSDAV to test 
case 1. Figure IV-23 presents frequency vs. power magnitude 
for three different time slice ensemble segments. All three 
segments were computed using the Hamming window and with the 
segment mean removed. The first ensemble segment, indicated 
by the solid line, is centered at sample 140, which is roughly 
the center of a portion of turbulent flow for the 41 pulses 
which comprise test case 1 (see Figure IV-1 for the relative 
position of the segments). The second ensemble segment, 
indicated by the short dashed line, is centered at sample 300, 
which is roughly the center of a portion of laminar flow. 
Observe the significantly lower power levels exhibited by the 
laminar segment at all frequencies above approximately 200 Hz. 
The third ensemble segment, indicated by the long dashed line, 
is centered at sample 500. Observe that, as expected, the 
power level has returned to that associated with turbulent 
flow. Figure IV-23 graphically illustrates the time dependent 


nature of the frequency spectra of non-stationary data. 
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The option to remove the mean was selected for 
this example, as with the stationary example discussed above, 
because the resulting zero frequency power spectrum value 
would otherwise be extremely large (on the order of 10) 
compared to the remainder of the frequency spectrum. The 
interesting part of the frequency spectrum is well above the 
propeller pulse frequencies where the turbulent fluctuations 
lie. Therefore the spectrum data below 200 Hz was not shown 
and the vertical scale was enlarged to better illustrate the 
high frequency portion of the spectrum. As discussed for the 
stationary case above, an alternate approach would be to use 
the low pass filtering algorithm (frequency domain smoothing) 
to compute and then remove the low frequency velocity 
component from each pulse in the ensemble (each column in 
array ENSMBL). This would be equivalent to high pass 
filtering. Then proceed with the algorithm without mean 
removal (since the "mean" has already been removed) to 
estimate the spectrum of the remaining high frequency 
fluctuating components of the ensemble segment. 

Examination of Figure IV-1, IV-2 and IV-3 suggests 
that the spectral characteristics of the unsteady boundary 
layer transition abruptly from laminar to turbulent flow and 
back again aS opposed to a gradual transition. If the 
ensemble segment crosses the transition point, the resulting 


spectrum is an average of the frequency content of the two 


106 


regions on either side of the transition point. A series of 
256 point wide ensemble segments moving across a transition 
point will produce a series of spectra which shows a gradual 
change from low power (laminar flow) to high power (turbulent 
flow). This phenomenon is due to the assumption of 
stationarity over a short time segment. If the segment 
encompasses a transition point, then the short segment 
stationarity assumption is violated (as discussed in paragraph 
D.4.a. above). The user should exercise care to insure that 
the ensemble segment does not cross a transition point. As 
a minumum, positioning the transition point near one edge of 
the ensemble segment will allow the tapering window function 


to minimize the effect of the abrupt transition. 


E. CORRELATION 


1. Overview 

The autocorrelation and cross correlation functions 
of stationary data were explicitly defined in Chapter III by 
equations III-17, III-20 and III-21. Equation III-32 presents 
an alternate definition of correlation (auto or cross) which 
is obtained indirectly through the use of the DFT. Algorithms 
are developed and implemented in this section for both the 
explicit and indirect methods. Then algorithms are discussed 
for the non-stationary case. Time limitations precluded 


implementation of the algorithms for non-stationary data. 
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2. Correlation of Stationary Data--Explicit Method 


a. Description of Algorithm 


The algorithm for computing either autocorrelation 


or cross correlation is essentially a direct implementation 


of equations III-17, III-20 and III-21. The algorithm is 


accomplished through the following steps: 


Identify the input and output data file names and specify 
the option to compute autocorrelation or cross 
correlation. 


If autocorrelation is selected, the input file must 
consist of a sequential file containing the sampled 
values of the random variable to be analyzed (e.g., 
instantaneous velocity, u or v). If cross correlation 
is selected, the input file must consist of a sequential 
file containing alternating values of the two random 
variables to be analyzed (e.g., instantaneous velocities 
u and v). The first variable of each pair of samples is 
the variable which is lagged. (For example enter 
alternating samples of u and v beginning with u to 


compute R,,). 

Select the number of lags (the maximum7 ) for which to 
compute the correlation. (Enter as a fraction of the 
length of the input data file). 


Remove the mean of the stationary data to obtain the 
fluctuating component. 


Compute the unbiased correlation function (auto or cross 
depending on the option selected). This step uses 
subroutine CORMAR which was extracted from Marple [Ref. 
101 

Write the output file. 


b. Implementation and Examples 


Figure IV-24 presents test case 2 which iS an 


example of steady (stationary) turbulence. This particular 


example was generated by measuring the velocity of the air 
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flow behind a grid or screen. The mean velocity was 95 
feet/second and the sample rate was 1000 Hz. 

Test case 2 was used as the example input to 
program CORREX which is the implementation of the explicit 
method for autocorrelation and cross correlation. The source 
code and user's guide for program CORREX are presented in 
Appendix B. Program CORREX was implemented to compute the 
correlation functions for positive lag values only. The 
autocorrelation function of a real time series is symmetric 
about the origin. Figure IV-25 presents the results of 
computing the autocorrelation function of test case 2 using 
program CORREX. As expected with random data, the 
autocorrelation of the fluctuating component at zero lag is 
a large value (the mean square of the fluctuating component). 
The autocorrelation function falls rapidly to oscillate about 
zero. In turbulence analysis, the autocorrelation function 
is integrated over all lags to obtain a measure of the 
"turbulence length scale" outlined in Chapter II. This 
essentially becomes the area under the initial peak. 

There is no actual two-component velocity data 
available with which to demonstrate the cross correlation 
option of program CORREX. Therefore two-component data were 
Simulated by using the single velocity component of test case 
2 twice, with the second velocity "component" (v) a time 


shifted version of the wu component. The "vy component" 
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Figure IV-24. TEST CASE 2 (AN EXAMPLE OF SsTATIONARY 
TURBULENCE. SAMPLE RATE--1000 HZ, MEAN VELOCITY--95.4 
FT/SEC.) 
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Figure IV-25. RESULT OF PROGRAM CORREX ON TEST CASE 2 
(AUTOCORRELATION OF STATIONARY TURBULENCE). 
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consists of the u component shifted by 50 samples. The 
expected result is a time shifted version of the 
autocorrelation with the peak occurring at a lag of 50 instead 
of zero. Illustrating the expected result, Figure IV-26 
presents the results of the cross correlation function of the 
simulated two component velocity signals computed using 
program CORREX. To obtain the cross correlation function at 
negative lag numbers, simply reverse the u and v input data 
vectors. In general, the cross correlation function is not 
symmetric. 


3. Correlation of Stationary Data--Fourier Transform 
Method 


a. Description of Algorithm 


The indirect or Fourier transform method for 
computing auto and cross correlation was outlined in section 
III.B.4. The idea is to compute the circular correlation per 
equation III-31, then normalize the circular correlation per 
equation III-32 to obtain the linear correlation. The input 
time series (u for autocorrelation, u and v for cross 
correlation) is zero padded prior to computing the circular 
correlation so that the correlation function at positive and 
negative lags may be separated. The algorithm is accomplished 
using the following steps: 

- Identify the input and output data file names and specify 


the option to compute autocorrelation or cross 
correlation. 
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Figure IV-26. RESULT OF PROGRAM CORREX ON TEST CASE 2 
(SIMULATED CROSS CORRELATION OF STATIONARY TURBULENCE). 
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If autocorrelation is selected, the input file must 
consist of a sequential file containing the sampled 
values of the random variable to be analyzed (e.g., 
instantaneous velocity, u or v). If cross correlation 
is selected, the input file must consist of a sequential 
file containing alternating values of the two random 
variables to be analyzed (e.g., instantaneous velocities 
uandv). The first variable of each pair of samples is 
the variable which is’ lagged. (For example enter 
alternating samples of u and v beginning with u to 


compute R,). 


Select the number of lags (the maximum 7 ) for which to 
compute the correlation. (Enter as a fraction of the 
length of the input data file). 


Remove the mean of the stationary data to obtain the 
fluctuating component. 


zero pad the data with at least as many zeros as the 
number of lags for which the correlation function will 
be computed (i.e., the number of zeros equals the maximum 
number of lags). If cross correlation was selected, then 
zero pad both input time series (e.g., u and Vv). 


Compute the circular correlation function (auto or cross 
depending on the option selected). This step uses 
subroutine CORREL which was extracted from Press, et al. 
(Ref. 12]. The circular correlation subroutine calls FFT 
routines which were also extracted from Reference 12. 
Note that this subroutine computes N times the result 
given by equation III-31 (where N is the number of points 
in the time series). Also note that both positive and 
negative lags are evaluated. 


Normalize the circular correlation by dividing the result 
of the previous step by the factor (N-r) where r is the 
lag number. This yields the desired linear correlation. 


Write the output file. 


b. Implementation and Examples 


The algorithm for computing the correlation 


functions using the indirect Fourier transform method was 


implemented as program CORFFT. The source code and user's 


guide are presented in Appendix B. Test case 2 was used as 


the example input to program CORFFT. The results are 
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virtually identical to those obtained from program CORREX, 
(see Figures IV-25 and IV-26) so no additional figures are 
presented. 

The primary advantage to program CORFFT is 
increased computational efficiency over the explicit method. 
The explicit method requires on the order of N° floating point 
operations to compute the correlation function. Using a FFT 
algorithm, the Fourier transform method computes’ the 
correlation function in roughly N log(N) floating point 
operations. As an example, using the explicit method (program 
CORREX), the autocorrelation example of test case 2 required 
26.7 seconds of computation time. A 50% time savings was 
achieved using the Fourier transform method )program CORFFT), 
which required only 13.4 seconds to obtain the same result. 
Even larger percentage time savings may be expected with 
longer input data files. (Test case 2 is only 1000 data 
points long.) Therefore program CORFFT is the recommended 
choice for computing correlation functions of stationary 
turbulence data. 

4. Correlation of Non-Stationary Data 
a. Description of Algorithm 

This section outlines some ideas on potential 
algorithms for computing correlation functions of non- 
stationary data. Time limitations prevented completion of the 


effort to develop these ideas into working computer programs. 
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Recall that the correlation function and the 
spectral density function form Fourier transform pairs. 
Therefore the algorithm for estimating non-stationary 
correlation functions should parallel that developed for non- 
stationary spectral estimation. One approach is proposed for 
computing autocorrelation. Simply compute the inverse DFT of 
the segment spectrum computed using program ENPSDAVG presented 


in section, 1V.C.4. 


However, a general algorithm for either 
autocorrelation or cross correlation is desired. Without 
repeating section IV.C.4., the basic idea is to divide the 


ensemble of non-stationary data into short time segments where 
the assumption of local stationarity may be applied. Positive 
lags up to half the length of the segment may be computed. 
The correlation function of the corresponding segment from 
each pulse is computed, then ensemble averaged to obtain the 
final correlation function for the ensemble of segments. 

As with the non-stationary spectral analysis, an 
important issue to consider in non-stationary correlation is 
the length of the segment. If the segment is too short, then 
a sufficient number of lags cannot be computed in order to 
characterize the turbulence length scale (autocorrelation). 
If the segments are too long, then the assumption of local 
stationarity is violated. At this point it is not clear that 
128 positive lags (corresponding to a 256 point wide segment) 


are sufficient to derive an accurate turbulence length scale. 
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or increasing the number of pulses so that a pure ensemble 
averaging algorithm may be applied. Again, the tradeoff is 
computer memory and storage limitations. 

Computer memory limitations become critical when 
cross correlation functions are desired. The array containing 
the ensemble of non-stationary data must now be twice as large 
to contain two input data vectors (u and v) instead of one. 
Rather than storing all of the data in computer memory at 
once, techniques may need to be developed to handle the data 
sequentially or in smaller blocks. 

The basic algorithm for computing non-stationary 
correlation functions is outlined in the following steps: 


¢ Enter the input data file name, the output data file name 
and the sample rate of the data. 


¢ Separate the input time history record into an ensemble 
of individual pulses. As previously described, each 
column of array ENSMBL represents a single propeller 
pulse. 


- Allow selection of one of three options for segmenting 
the ensemble. The first option segments the ensemble 
based on an input value for the number of data points 
between the beginning of each segment (i.e., the amount 
of time between each PSD). The second option segments 
the data based on an input value for the number of 
segments across the ensemble of pulses (i.e., the number 
of PSDs across the ensemble). The third option creates 
a Single segment ensemble based on an input time (sample) 
value of the center of the segment (i.e., choose the 
location for a single PSD). 


* Remove the mean of each segment. Then zero pad each 
segment such that the positive and negative lags may be 
separated. 


¢ Compute the circular correlation estimate of the 
corresponding segment of each pulse. Average the 
correlation function from each pulse segment and 
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normalize to obtain linear correlation from circular 
correlation. 


+ Write the output file. 
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V. CONCLUSIONS AND RECOMMENDATIONS 


A data analysis system was developed to support a planned 
wind tunnel investigation into the nature of unsteady 
turbulent boundary layers. A series of algorithms were 
developed and implemented to compute the statistical 
parameters which the aerodynamicist uses to characterize the 
nature of turbulent fluid § flow. Specifically, the 
aerodynamicist will use these data analysis systems to analyze 
fluid velocity measurements obtained in the boundary layer of 
an airfoil subjected to periodic turbulent wake disturbances 
such as in a propeller slipstream. The statistical parameters 
used to characterize the unsteady turbulent boundary layer 
include the mean velocity, turbulence intensity, power 
spectral density function and the autocorrelation or cross 
correlation functions. Unsteady turbulent fluid flow is a 
non-stationary random process. Therefore the statistical 
parameters of interest are time varying functions and require 
special non-stationary techniques to estimate. 

The following paragraphs summarize the specific algorithms 
included in the data analysis system. Where alternate 
algorithms are proposed for a single statistical parameter, 


the "best" method is recommended base on specified criteria. 
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A. MEAN VELOCITY 


Four different algorithms were developed and implemented 
to determine the non-stationary mean velocity at a point in 
the boundary layer. Program ENSEMBL estimates the non- 
stationary mean by performing ensemble averaging on an 
ensemble of propeller blade wake passages. Storage and 
computer memory limitations preclude collecting enough data 
to obtain an adequate estimate of the mean using only ensemble 
averaging. Program MOVAVE performs a moving time average on 
each pulse in the ensemble to obtain a "local mean" prior to 
ensemble averaging. Program SMUMEAN performs frequency domain 
low pass filtering on each pulse to obtain a "local mean" 
prior to ensemble averaging. Program MOVAVE and SMUMEAN both 
yield improved estimates of the non-stationary mean. Program 
SMUMEAN achieves the following advantages: 


+ Increased flexibility due to a large number of available 
cutoff frequencies. 


- Ability to implement alternate low pass filters, if 
desired. 


- No loss of data at the beginning and end of the ensemble 
average data record. 


Program MEANSMU performs the ensemble averaging operation 
prior to frequency domain low pass filtering. This achieves 
the identical result as program SMUMEAN with an approximately 
75% reduction in computation Sie. Therefore, program MEANSMU 


is the recommended method for determining the non-stationary 
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mean velocity. Program MEANSMU (and SMUMEAN) was implemented 
with only one type of low pass filter. Therefore, it is 
further recommended that alternate frequency domain low pass 
filters be implemented and tested in the smoothing algorithm 
of program MEANSMU to determine if there is any significant 
difference or improvement in the non-stationary mean velocity 
estimate. In particular, filters which better approximate the 


ideal low pass filter should be examined. 


B. TURBULENCE INTENSITY 


Turbulence intensity 1S a measure of the amount of 
variation of the fluctuating velocity component about the 
local mean velocity. Two different algorithms were 
implemented based on two methods for estimating the local 
mean. Program TURINTMA computes the local mean via a moving 
average. Program TURBIN computes the local mean via frequency 
domain low pass filtering. Programs TURINTMA and TURBIN 
produce nearly identical turbulence intensity results because 
the frequency response characteristics of the two different 
smoothing algorithms are similar. The smoothing algorithms 
are the same as those used in the mean velocity algorithms, 
so the same advantages described above apply to program 
TURBIN. Therefore, program TURBIN is the recommended method 
for estimating non-stationary turbulence intensity. It is 


further recommended that alternate frequency domain low pass 
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filters be implemented and tested in the smoothing algorithm 
of program TURBIN to determine if there 1s any significant 
difference or improvement in the non-stationary turbulence 
intensity estimates. In particular, filters which better 


approximate the ideal low pass filter should be examined. 


C. SPECTRAL ANALYSIS 


Two programs were developed to estimate the spectral 
characteristics of the turbulent boundary layer velocity data. 
Program PSD was developed for steady (stationary) turbulence. 
The stationary case waS examined first to gain an 
understanding of the techniques involved. Program PSD is an 
implementation of the classical Welch "overlapping segments" 
periodogram method of spectral estimation. The Welch method 
1S non-parametric and as such, does not depend on the 
assumption of a model for the random data in order to obtain 
accurate spectral estimates. Several spectral window function 
options are included for leakage reduction. Window option 5 
(Hamming window) or 6 (3 term Blackman-Harris window) are the 
recommended choices of window functions for the generally 
broad band turbulence data. These window functions provide 
a good compromise between the ideals of narrow. main lobe width 
and low side lobe level. However, due to the broad band 
nature of the data, any of the popular windows (options 4 


through 7 in Table IV-4) will produce similar results. 
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Program PSD also provides an option for removal of the 
stationary mean prior to estimating the spectrum. This option 
should be selected for data with a large mean value. Due to 
Spectral leakage (which windowing reduces but does not 
eliminate) the large zero frequency power may mask the true 
power of nearby (low) frequencies. It is further recommended 
that an alternate option for mean removal be implemented. 
Rather than removing only the mean value, an option should be 
implemented which will high pass filter the data and remove 
all low frequency content below the frequencies of interest. 

Program ENSPSDAV was developed to estimate the time 
varying non-stationary spectral characteristics of an ensemble 
of propeller blade wake passages. This program computes a 
spectral estimate for a short time segment of each propeller 
pulse in the ensemble. The data are reasonably assumed to be 
stationary over a short segment. The program then ensemble 
averages the individual spectral estimates to obtain the final 
estimate for the short time segment of the ensemble. The non- 
stationary (time varying) spectral characteristics are 
constructed by examining different short time segments. 

The same options for window functions and mean removal are 
provided in program ENSPSDAV as are provided in program PSD 
discussed above. The same recommendations for window function 
choice and an alternate mean removal method also apply to 


program ENSPSDAV. It is further recommended that the user of 
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program ENSPSDAV exercise care in the selection of the short 
time segment locations. The user should insure that the short 
time segment does not cross a transition point between laminar 
and turbulent flow, thus violating the assumption of 
stationarity over the short segment. If a short segment 
cannot be located without including a transition point, then 
the transition point should be located near the beginning or 
end of the segment such that the window function will minimize 


the effect of the abrupt transition point. 


D. CORRELATION 


Two programs were developed to estimate the 
autocorrelation or cross correlation of stationary data. 
Program CORREX is a direct implementation based on the 
explicit definition of autocorrelation and cross correlation. 
Program CORFFT estimates the autocorrelation or cross 
correlation using the indirect Fourier transform method. 
These two programs produce virtually identical results. 
Program CORFFT has the primary advantage of increased 
computational efficiency. For a test case data record which 
is 1024 points in length, a 50% computational time savings was 
achieved. Greater percentage savings are expected for longer 
data records. Therefore, program CORFFT is the recommended 
choice for estimating correlation functions of stationary 


turbulence data. 
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Some ideas were presented for potential algorithms to 
estimate the correlation functions of non-stationary 
(unsteady) turbulence data. Time limitations prevented the 
development of these ideas into working computer programs. 
It is recommended that the algorithms for estimating the non- 
stationary autocorrelation or cross correlation functions 
presented in Chapter IV be developed into working computer 


programs. 
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APPENDIX A 


FREQUENCY RESPONSE CHARACTERISTICS OF 
THE SMOOTHING ALGORITHMS 


A. FREQUENCY RESPONSE OF THE MOVING AVERAGE SMOOTHING 
FILTER 


A moving average smoothing filter is one form of a non- 
recursive time domain low pass digital filter. Moving average 
smoothing involves computing a time average over a short 
segment of a time series. The average of the short segment 
represents the smoothed value for the center point of the 
segment. This is equivalent to the midpoint of the linear 
least squares curve fit through the segment of points. For 
a discrete time series u(t,) with equally spaced samples, the 


moving average filter is given by the equation: 





N 

= =e A-1 

where u. is the smoothed value or "local mean" 
1LOCAL 


corresponding to the instantaneous value u(t,). L is the 
length of the moving average window (an odd number) such that 
L= 2N+ 1. In other words, N is the number of points on each 
side of the i'th "pivot point". The index m runs from -N to 


Ne Retmon 
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To examine this digital filter in the frequency domain, 
examine the transfer function of the digital filter. The 
transfer function of the digital filter is the eigenvalue of 
equation A-1. To find the transfer function (eigenvalue), 
substitute an eigenfunction for the input time series u(t,), 
then solve for the eigenvalue. The trigonometric sinusoids 
are the eigenfunctions of an equally spaced sampled process 
and form an orthogonal basis spanning the finite time series 
interval [Ref. 9]. Writing the trigonometric sinusoids in 


complex exponential form, let 


Mie =<! (A-2) 


where 4) is the frequency in radians. Substituting the 


eigenfunction into equation A-1 yields: 
H(w) = 55 S elmo (A-3) 
2N+1 Sy 


where H ( &) is the transfer function. Note that equation 
A-3 looks just like the discrete Fourier transform with an 


input function of 1. Expand the summation: 


H(w) = seer eNO I, IN 2s tt. telh Op eNO) (a4) 


But the exponential terms may be expressed as sines and 
cosines as follows: 
a! = Cos(-Nw) +jsin(-Nw) 


e M00 = cos/-(N-1) 0} + jsin[-(N-1)o] 
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et 0 = cos(-(N-2) a] =h jsin[-(N-2)0] 


JN Do = cos| (N-1)a] ~ jsin| (N-1)o) 


N oe 
e/Xo = cosNw + jsinNw 


Substituting these identities into equation A-4 and 
Simplifying yields the transfer function for the moving 


average filter: 


sin(N+5) w 
H(w) = eee ee G) (A-5) 
(2N+1) sins 


Using equation A-5 and the identity: 


w = an{ =| (Asay 


where f is the frequency in Hertz and oe is the sample rate, 
the frequency response characteristics of the moving average 
filter may be computed and plotted for different window sizes 
(L = 2N + 1). Figure A-1 presents the transfer function of 
the moving average filter for several different window sizes 
at frequencies up to the Nyquist frequency (see section 
III.B.2.). The cutoff frequency of each filter is defined as 
the frequency anore the magnitude has decreased to 0.7071 (3 
dB) of its initial value. Table IV-1 presents the cutoff 
frequencies for many different window sizes. These were 
determined by solving the following equation for f,, the 
cutoff frequency: 
sin| (w+3) 20( #4) 
S 


0.7071 = (A-7) 
(2N+1) sina | ¢*] 
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Figure A-1. 
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Program FCMA is provided in Appendix B to compute the solution 
for different sample rates and window sizes. 

Observe the large oscillations or ripples in the frequency 
response characteristics of the moving average filter. These 
ripples are known as the Gibbs' phenomenon and are caused by 
the truncated Fourier series of equation A-3 [Ref. 9]. 
Truncating the Fourier series is necessary for a finite time 
series. The Gibbs' phenomenon causes the filter transfer 
function to yield a poor approximation of the ideal low pass 
filter presented in Figure IV-9. Finding better 
approximations to the ideal low pass filter provides the 


motivation for examining other smoothing algorithms. 


B. FREQUENCY RESPONSE OF THE FREQUENCY DOMAIN SMOOTHING 
FILTER 


Instead of operating on the time series in the time 
domain, a more natural method to accomplish low pass filtering 
(smoothing) is to transform the time series into the frequency 
domain using the DFT. Then the input time series (now a 
"frequency series") is simply multiplied by a low pass filter 
whose transfer function is defined in the frequency domain to 
possess the desired frequency response filter characteristics. 
Transfer functions may be defined which minimize or eliminate 
the Gibbs' phenomenon ripple effect and better approximate the 


ideal low pass filter. 
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One such transfer function, used in the frequency domain 
smoothing filter described in Chapter IV, is a parabolically 


shaped filter defined by: 


H( +] = wax |1-(f-7) ,0| J=0,1,2,...,M (A-8) 
Ss 


where L is the "smoothing window" size which determines the 
filter cutoff frequency, M is the (integer power of two) 
number of data points in the time series, and J is the 
frequency bin index running from 0 to M. Figure A-2 
summarizes the frequency response characteristics of this low 
pass filter. The plot presents J versus magnitude for several 
different smoothing window sizes. The window size is 
parameterized as L/M since the actual frequency response 
characteristics are a function of the number of points in the 
time series. The actual frequency associated with the index 


J 1S given by: 
f=f i) 
s| M (A-9) 
As before, the cutoff frequency is defined as the 
frequency at which the magnitude has attenuated to 0.7071 of 


the initial value. Combining equations A-7 through A-8, the 


cutoff frequency may be determined by solving the equation: 


0.7071 = 1-[u( #4 | (A-10) 
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Which yields: 


f. = o.2929( 42) (A-11) 


Table IV-2 was generated using equation A-11. 


C. COMPARISON OF THE TRANSFER FUNCTIONS 


Figure A-3 presents a comparison of the transfer functions 
of the moving average smoothing filter and the frequency 
domain smoothing filter. Figure A-3 is presented for a sample 
rate of 15,000 Hz and a time series length of 1024 points. 
The smoothing window widths were selected such that the cutoff 
frequencies of the two filters are nearly equal. For the 
moving average filter, f, = 320 Hz, corresponding to a window 
size of 21 points. For the frequency domain filter, f, = 325 
Hz, corresponding to a 25 point smoothing window. The two 
transfer functions are nearly equal at frequencies up to the 
cutoff frequency. The frequency domain filter falls off more 
rapidly and does not oscillate. Therefore, the frequency 


domain smoothing filter provides a better approximation to the 


ideal low pass filter. 
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APPENDIX B 


SOURCE CODE AND USER’S GUIDES FOR FORTRAN 
PROGRAMS 


A. GENERAL INFORMATION 

This section provides general information and guidance 
which applies to all of the programs presented in this 
Appendix. Subsequent sections present user documentation 
specific to each individual program. 

The data acquisition system planned for the tests 
described in Chapter II is an IBM PC/AT (R) clone 
microcomputer equipped with an analog-to-digital conversion 
board and an 80287 math co-processor. The data will be 
Gigitized and recorded on the hard disk in real time. Data 
analysis using the programs described in this thesis will be 
performed subsequent to the tests. These programs are 
intended to run on the same hardware setup just described. 

All of the programs presented in this thesis were written 
and compiled with Microsoft (R) FORTRAN version 4.1 operating 
under Microsoft's MS DOS (R) version 3.2 operating system. 
Microsoft's (R) FORTRAN compiler implements the full ANSI 77 
standard for the FORTRAN programming language. This version 
of FORTRAN also includes many extensions to ANSI 77 standard 


FORTRAN. Compiler options were selected to take full 
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advantage of the 80287 math co-processor. All program and 
data file names follow the MS DOS (R) standard format for file 
names. 

The programs presented in this Appendix are all heavily 
commented and documented to allow the user to follow the logic 
and make modifications as necessary. Other individual 
applications will likely have different input file formats. 
So, the input file "read" statements are most likely to need 
modification. Most of the programs presented here have 
integrated the input file "read" statements into the algorithm 
which separates the input velocity data into an ensemble of 
pulses. This algorithm will be described first and applies 
to all programs which operate on the ensemble of pulses. 
(Programs PSD, CORREX, CORFFT and MEMPSD do not operate on an 
ensemble.) 

The input file is expected to be an ASCII sequential file 
with alternating values of two parameters (either in a single 
column or two adjacent columns). The two parameters are 
expected to be velocity (U) and the propeller trigger (PROP), 
in that order. The propeller trigger signal provides a sharp 
spike at each propeller blade passage and is used to separated 
the velocity data into an ensemble. The program reads the 
velocity and prop trigger data until the first trigger spike 


is encountered. This is the read statement that should be 
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modified to handle other input file formats. The first 
trigger spike exceeding a value of one (volt) initiates the 
beginning of the first blade wake passage. The velocity 
samples are copied into the first column of array ENSMBL as 
they are read from the input file. Upon encountering the next 
trigger spike, the velocity samples are copied into the 
beginning of the next column. In this manner, array ENSMBL 
is filled with velocity samples such that each column 
represents a blade wake passage. Data is read until 50 blade 
wake passages have been copied into array ENSMBL or until an 
end-of-file (EOF) is encountered. (The last partial pulse 
which was being read when the EOF was encountered is deleted 
since it is incomplete.) The length of the shortest column 
is saved (LASTI) so that the entire array may be truncated to 
the length of the shortest blade wake passage. 

The length of each column is displayed on the screen as 
the program runs. The user should watch these values 
carefully to insure that the columns are all about the same 
length. Any significant difference in pulse length (more than 
about five to ten samples) is a possible indication of a 
problem with the propeller trigger signal. Erroneous results 
may be obtained if there are significant discrepancies in the 


length of the blade wake passages. 
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B. 


USER’S GUIDE FOR PROGRAM ENSEMBL 


1. Input Instructions 


Insure that the executable file for program ENSEMBL 
(ENSEMBL.EXE) and the input data file are in the same 
subdirectory of a disk drive with sufficient space 
remaining to write the output file. 


Initiate program ENSEMBL (type "ENSEMBL" at the DOS 
prompt). 


The program will prompt the user for the name of the 
input data file. Enter the full name of the input data 
file including the extension. 


The program will prompt the user for the name of the 
output file which will be created by the program. Enter 
the desired file name. Do not use the name of an 
existing file or the program will terminate. 


The program will then read the input data file, compute 
the ensemble average and write the output data file. 
Remember to observe the column lengths of the array 
ENSMBL to insure that the ensemble is being constructed 
properly. 


2. Output Description 


Program ENSEMBL computes the ensemble average and 


writes an output file consisting of three columns. The first 


column is ensemble sample number (or row number, running from 


1 to LASTI). The second column is the ensemble average 


values. The third column is the ensemble variance. The 


following pages present the source code for program ENSEMBL. 


138 


PROGRAM ENSEMBL 


KREREEKAKAARKAKEKARKKREKEAAEKAARKEEREEEKAKKKKAKKKKRKKKKRKRKKRKRKRKEKEKKEKKRKKKRKEKKSE 


i A A 


Program to perform ensemble averaging on an ensemble of unsteady 
turbulent boundary layer velocity measurements such as occurs 
in the boundary layer of an airfoil immersed in a propeller slipstream 


Input data file is assumed to be in output format of the A/D conversion 
system. For this application the input file is a sequential 

file with alternating values of velocity (U) and propeller pulse 
(PROP), starting with velocity. 


The program first separates the single time history input file into an 
ensemble of pulses using the value of the propeller trigger (PROP) 

as the marker for the beginning of a pulse. The program also calculates 
statistical variance and a "pseudo" turbulence intensity for the ensemble. 


Key variables are defined below: 


ENSMBL = The primary data array which contains up to 50 pulses (columns) 
each of length up to 2048 points. The pulses are synchronized in 
time (an ensemble) such that the beginning of each column (i.e. 
the first row) represents the beginning of each pulse. 

U = The measured velocity as read from the input file into array ENSMBL. 

PROP = The propeller pulse timing marker signal. Also read from the input. 

IPLS = Number of pulses to be read from the input file. 

LASTI = Number of points in the columns of ENSMBL. (ie. rows in ENSMBL) 

NPLS = (J-1) = Final total number of pulses (columns in ENSMBL). 

AVG = The ensemble average of U. (Output) 

VAR = The ensemble statistical variance of U. (Output) 

TURBIN = The "pseudo" ensemble turbulence intensity of U. (Not based 
on local smoothed velocity but on the overall ensemble mean 
velocity. Really a variance.) Output. 

UINF = Free stream velocity external to the boundary layer. 


Programmed by: Donald K. Johnson 
July 1988,NPS 


RRA ERRERKRRAKKRAKKRRERRRRREKERRRAERERAEARREKARERKAKRKKKKEKRKAKKKEKKRKKKKKKRKKKK 


DIMENSION ENSMBL(2048,50) ,AVG(2048) , VAR(2048) , TURBIN (2048) 
INTEGER I,IPLS, LASTI,J,K,L,NPLS 
CHARACTER*12 INFILE, OUTFILE 


PRINT'(/A/)',' Program to compute pure ensemble average:' 


PRINT '(/A\)',' ENTER THE INPUT DATA FILE NAME: ' 

READ (*,'(A)') INFILE 

PRINT '(/A\)',' ENTER THE NUMBER OF PROP PULSES YOU WANT: ' 
READ (*,*) IPLS 

IPLS=50 


PRINT '(/A\)',' ENTER THE OUTPUT DATA FILE NAME: ' 
READ (*,'(A)') OUTFILE 


£39 


* ee * * 


* 


x 


10 
akt 


kkk 


zak 


20 


28 


1mrM. 


pseudo turbulence intensity calcs are currently de-activated *## 
PRINT '(/A\)',' ENTER U-edge (REFERENCE VELOCITY FOR TURB.INTENSITY 
1 NORMALIZATION): ' 

READ(*,*) UINF 

UINF=UINF*1. 

UINF=1. 


OPEN(10, FILE=INFILE, STATUS='OLD') 
OPEN (11, FILE=OUTFILE, STATUS='NEW') 


LASTI=2048 

J=1 

I=0 

PRINT'(/A/)",' READING INPUT FILE....° 

PRINT'(/A/)',' WARNING, Watch the number of points in each column 


1... they should all be about the same length.' 


CONTINUE 


Read the data into array "ENSMBL", each column is a pulse *** 
READ (10,*,END=20,ERR=20) U, PROP 
I=I+1 
But skip the first partial pulse ****s* 
IF (J .GT. 1) ENSMBL(I,J<-1) =U 
If you have a new pulse, start a new column: 
a new pulse is defined by PROP being greater than 1. volt 
IF(I .GE. 20 .AND. PROP .GT. 1.) THEN 

MF (9° .EO. 1) THEN 

WRITE(*,*) 'Skipped first',I,' points (not a full pulse).' 
ELISE 
WRITE(*,*) ‘No. Of points in column’; J=-1, sis ,1 

ENDIF 

Find the shortest column *** 

IF(I .LT. LASTI .AND. J .Gioeteieasit.—. 

Quit loop if you have enough pulses ***# 

IF (J-1 .GE. IPLS) GO TO 28 

J=J+1 

I=0 
ENDIF 


GO TO 10 
CONTINUE 


WRITE(*,*) 'End of File encountered - deleting last partial colu 
t 


J=Jol1 


CONTINUE 
XN=FLOAT (J-1) 
NPLS=J-1 


WRITE(*,*)'The final number of pulses is',NPLS,' each of length', 
1 LASTI 


zzz Compute the ensemble mean, variance and "pseudo" turbulence intensity 


* 
Cc 


based on user input value of V - infinity (free stream velocity). 


DO 22 I=1,LASTI 


SUM=0. 
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SUMSQ=0. 
DO 30 L=1,NPLS 
ba lalial first compute the sum *** 
SUM2SUM+ENSMBL(I,L) 
SUMS Q=SUMSQ+ENSMBL(I,L) **2 
30 CONTINUE 
alia then divide by the number of points 
AVG (I) =SUM/XN 
VAR (I) = (SUMSQ~-XN*AVG (I) **2)/ (XN—-1.) 
TURBIN (I) =(( (SUMSQ-XN*AVG (I) **2) /XN) **.5) /UINF 
22 CONTINUE 
x*e* Write the output *****e% 
DO 40 I=1,LASTI 
XI=FLOAT (I) 
Eel note: the pseudo turbulence intensity is currently not output 
WRITE(11,200) XI,AVG(I),VAR(I) 
200 FORMAT(F5.0,F8.1,F6.1) 
40 CONTINUE 
CLOSE (10) 
CLOSE (11) 
END 
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C. USER’S GUIDE FOR PROGRAM MOVAVE 


1. Input Instructions 


* Insure that the executable file for program MOVAVE 
(MOVAVE.EXE) and the inpute data file are in the same 
subdirectory of a disk drive with sufficient space 
remaining to write the output file. 


* Initiate program MOVAVE (type "MOVAVE" at the DOS 
prompt). 


*- The program will prompt the user for the name of the 
input data file. Enter the full name of the input data 
file including the extension. 


* The program will prompt the user for the name of the 
output file which will be created by the program. Enter 
the desired file name. Do not use the name of an 
existing file or the program will terminate. 


- The program will prompt the user for the moving average 
window size. The window size specifies the number of 
samples which are time averaged to produce an interim 
"local" mean or "local" smoothed velocity which is then 
ensemble averaged. The window size should be an odd 
number. Table IV-1 in Chapter IV presents some examples 
of the relationship between moving average window size 
and its associated low pass filter cutoff frequency. 
Other cutoff frequencies may be determined by solving 
equation A-7 in Appendix A. Program FCMA is provided 
(immediately after the source code for program MOVAVE) 
to solve equation A-7. 


- The program will then read the input data file, compute 
the non-stationary mean and write the output data file. 
Remember to observe the column lengths of the array 
ensembl to insure that the ensemble is being constructed 
properly. 

2. Output Description 
Program MOVAVE computes the non-stationary mean and 


writes an output file consisting of two columns. The first 


column is ensemble sample number (or row number, running from 
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1 to LASTI). Knowing the sample rate, the time scale may be 
inferred from sample number. The second column contains the 
non-stationary mean values. The following pages present the 


source code for program ENSEMBL. 


a 


SDEBUG 


PROGRAM MA 
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Program to perform simple moving average smoothing on an ensemble of 
propeller pulses then ensemble average the smoothed data. Note that 

the moving average operation sets the first (L-1)/2 points and the last 
(L-1)/2 points of the mean velocity output to zero (where L is the 
smoothing window size). The relatioship between smoothing window 

size and cutoff frequency is given by: 


SIN(((L-1)/2 + .5)*2*PI*FC/FS) / L*SIN(PI*FC/FS) = 0.7071 


where FC is the cutoff frequency, FS is the sample rate and L is the 
odd numbered window size. The table below gives some examples: 


ie FC (FS=30000 Hz) 
3° 4660 Hz 
11 1220 
17 790 
21 640 
35 380 


Input data file is assumed to be in output format of the D/A 
conversion system, il.e. a sequential file with alternating values 
of velocity (U) and propeller pulse (PROP), starting with velocity. 
The program first separates the single time history file into an 
ensemble of pulses using the value of the propeller trigger (PROP) 
as the marker for the beginning of a pulse. 


Key variables are defined below: 


ENSMBL = The primary data array which contains up to 50 pulses (columns) 
each of length up to 2048 points. The pulses are synchronized in 
time (an ensemble) such that the beginning of each column (i.e. 
the first row) represents the beginning of each pulse. 

U = The measured velocity as read from the input file into array ENSMBL. 

PROP = The propeller pulse timing marker signal. Also read from the input. 

IPLS = Number of pulses to be read from the input file. 

LASTI = Number of points in the columns of ENSMBL. (ie. rows in ENSMBL) 

NPLS = (J-1) = Final total number of pulses (columns in ENSMBL). 

L = The width of the moving average window (should be an odd number). 

A(I) = The current column of ENSMBL to be smoothed using subroutine MOVAVE 

B(I) = The current smoothed column of ENSMBL. 

USMTH(I) = The ensemble average of all the smoothed propeller pulses. 


Programmed by: Donald K. Johnson 
AUGUST, 1988,NPS 


kekkkekakekekkkekkkkkkkkkkkkkkkkkkkkkkkkkkkkkakkekakkeakkekkkaxkkkeakkkdkk&keekkkeaekkee 


PARAMETER (MMAX=2048) 
DIMENSION ENSMBL(2048,50) ,A(MMAX) , B(MMAX) , USMUTH (MMAX) 
INTEGER L,NPTS,ILAST,NPLS,IPLS . 

CHARACTER*12 INFILE, OUTFILE 
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10 
ake 
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PRINT'(/A/)',' Moving average smoothing and ensemble averaging. ' 


PRINT '(/A\)',' ENTER THE INPUT DATA FILE NAME: ' 
READ (*,'(A)') INFILE 


PRINT '(/A\)',' ENTER THE OUTPUT FILE NAME: ' 
READ (*,'(A)') OUTFILE 


PRINT '(/A\)',' ENTER THE MOVING AVERAGE WINDOW SIZE" (AN ODD NUMB 
1ER): ' 
READ (*,*) L 


OPEN (UNIT=10, FILE=INFILE, STATUS='OLD') 
OPEN (UNIT=11, FILE=OUTFILE, STATUS='NEW') 


separate the time history data into an ensemble of pulses **#* 
LAST I=2048 


IPLS=50 

J=1 

I=0 

PRINT'(/A/)',' READING INPUT FILE....! 

PRINT'(/A/)',' WARNING, Watch the number of points in each column 
1... they should all be about the same length.' 

CONTINUE 


Read the data into array "ENSMBL", each column is a pulse *** 
READ (10,*,END=20,ERR=20) U,PROP 
I=I+1 
But skip the first partial pulse *###** 
IF (J .GT. 1) ENSMBL(I,J-1) =U 
If you have a new pulse, start a new column. *** 
Note, prop trigger is set to 1. volt to start a new pulse **#% 
IF(I .GE. 20 .AND. PROP .GT. 1.) THEN 
TF (0 .EQ. 1) THEN 
WRITE(*,*)'Skipped first',I,' points (not a full pulse).' 
ELSE 
WRITE(*,*)'No. of points in column',J-1,' is’,I 
ENDIF 
Find the shortest column *** 
PCL. ie.s GASTE «AND. J .GT. 1) LAST#E#I 
Quit loop if you have enough pulses *** 
IF (J-1 .GE. IPLS) GO TO 28 
J=J+1 
I=0 
ENDIF 
GO TO 10 
CONTINUE 
WRITE(*,*)'End of File encountered - deleting last partial colu 
inte 
J=J-1 
CONTINUE 
XN=FLOAT (J-1) 
NPLS=J-1 
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WRITE(*,*) 'The final number of pulses is',NPLS,' each of length', 
1 LASTI 


apply the moving average window to each column (pulse) in array ENSMBL 
writing the averaged value back into array ENSMBL at the current point. 


PRINT '(/A/)',' WORKING .....' 
start the pulse loop *** 
DO 25 K=1,NPLS 
and start the time loop *** 
DO 15 I=1,LASTI 
copy the current (Kth) pulse into A for call to subroutine *** 
A(I)=ENSMBL(I,K) 
CONTINUE 
smooth the Kth pulse *** 
CALL MOVAVE(A, LASTI,L,B) 


write the smoothed values of velocity back into ENSMBL 
DO 16 I=1,LASTI 
ENSMBL (I,K) =B(I) 
CONTINUE 
CONTINUE 


now ensemble average (across rows) to get the final mean velocity 
start with a time loop *** 
DO 50 I=1,LASTI 
now sum across each row ***% 
SUMU=0. 
DO 60 K=1,NPLS 
SUMU=SUMU+ENSMBL(I,K) 
CONTINUE 
compute the ensemble average *** 
USMUTH (I) =SUMU/NPLS 
CONTINUE 
write the output *** 
DO 70 I=1,LASTI 
XI=FLOAT (TI) 
WRITE(11,100) XI,USMUTH(T) 
FORMAT (1X,F5.0,1X,F10.2) 
CONTINUE 
CLOSE(10) 
CLOSE (11) 
END 


CARRERE KKRKEKKRARKEKKEKRKEKKKKKKKKEKKEKKKKEKKKKKKE 
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SUBROUTINE MOVAVE(A,N,L,B) 


Program performs a simple moving average of a time series to determine a 
local mean (smoothed) velocity. Note that the first L/2 values and the 
last L/2 values of the input time series are set to zero. Briefly, the 
moving average accomplishes a crude low pass filter type of smoothing 

on the input time series as follows. Starting with the first L (usually 
odd) points in the time series, a time average is computed for the L 
points. This average value becomes the smoothed velocity at the center 

of the L points, point (L=-1)/2 for odd L. Then the moving average window 
is moved over one point and a new average is computed for the window. This 
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is repeated until the Lth point of the window reaches the end of the time 
series. 


The key variables are defined below: 


= A column array containing the time series to be averaged. 

= The number of values in A. 

= The size of the averaging window (an odd number is best because 
then the average of the small time segment is centered on the 
mid-point of the segment). Note, a larger window size corresponds to 
a lower low pass filter cutoff frequency. See the table below for the 
relation between window size and cutoff frequency. 

B = A column array containing the smoothed (averaged) time series values. 


oz 


DIMENSION A(1),B(1) 

AL=L 

MH=L/2 

zero out the answer array *** 
DO 18 I=1,N 


18 B(I)=0. 


K=L-1 
if L is .GT. N then quit *** 
cP SCN=15) 72) 2,2 
2 NN=N-K 
start the main averaging loop, working on L/2 points on either 
side of the KK th data point. (KK is calculated below) 
DO 4 J=1,NN 
JN=J+K 
reset the accumulator *** 
S=0. 
sum across L points centered on point KK. 
DO 3 I=J,JN 
S=S+A(I) 
3 CONTINUE 
now divide by L to get the average of segment centered on point KK. 
KK=J +MH 
B(KK)=S/AL 
4 CONTINUE 
if L is even then take the average of two adjacent terms to get the value 
of the average velocity centered on the KK th point. 
ieneo-(L/2)*2) 1,6, 1 : 
6 MS=MH+1 
MD=MH+NN-1 
DO 7 I=MS,MD 
TEMP=(B(I)+B(I+1))/2. 
B(I)=TEMP 
CONTINUE 
RETURN 
END 


rw) 
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PROGRAM FCMA 
Computes the relationship between the moving average smoothing window 
Size (L, where L=2*N + 1), the sample rate (FS), and the frequency (F), 
so that the low pass moving average smoothing filter cutoff frequency 
may be determined. Cutoff frequency is defined as the frequency where 
the magnitude is -3 dB (a factor of 0.7071) down from the unfiltered 
magnitude. 


PI=4.*ATAN(1.) 
PRINT '(/A\)',' ENTER SAMPLE RATE: ' 
READ(*,*) FS 


start the window size loop (window size, L = 2N+1l) **#* 
this loop will compute the cutoff frequency for all window sizes 
between 3 and 61. 
DO 10 N=1,30 
start with 5 Hz. frequency and step up in 5 Hz. steps *** 
F=5. 
KNT=0 
CONTINUE 
compute the magnitude *** 
XN=FLOAT (N) 
X1=SIN((2*PI*F/FS) * (XN+.5) ) 
X2=(2.*XN+1.) *SIN(PI*F/FS) 
H=X1/X2 
L=2*N+1 
if the magnitude is near the cutoff point then print results. *** 
the correct cutoff frequency is the one where the magnitude is 
less than or equal to 0.7071 
TFCH. ~LT. 0.708) THEN 
WRITE(*,100) L,F,H 
FORMAT (1X,1I5,F8.0,F10.6) 
KNT=KNT+1 
after writing results for 3 frequencies, skip to next window size 
EPORKNT: ~<GT. 3) “GOsTO. 410 


ENDIF 
increment frequency *** 
F=F+5. 
IF (F .LT. FS) GO TO 20 
CONTINUE 
END 
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D. USER’S GUIDE FOR PROGRAM SMUMEAN 


1. Input Instructions 


¢ Insure that the executable file for program SMUMEAN 
(SMUMEAN.EXE) and the input data file are in the same 
subdirectory of a disk drive with sufficient space 
remaining to write the output file. 


¢ Initiate program SMUMEAN (type "SMUMEAN" at the DOS 
prompt). 


¢ The program will prompt the user for the name of the 
input data file. Enter the full name of the input data 
file including the extension. 


¢- The program will prompt the user for the name of the 
output file which will be created by the program. Enter 
the desired file name. do not use the name of an 
existing file or the program will terminate. 


¢ The program will prompt the user for the FFT smoothing 
window siZe. The window size specifies the cutoff 
frequency of the low pass frequency domain filter. The 
low pass filtered data produces an interim "local" mean 
or local smoothed velocity which is then ensemble 
averaged. Table IV-2 in Chapter IV presents some 
examples of the relationship between moving average 
window size and its associated low pass filter cutoff 
frequency. Other cutoff frequencies may be calculated 
using equation A-11 in Appendix A. 


- The program will then read the input data file, compute 
the non-stationary mean and write the output data file. 
Remember to observe the column lengths of the array 
ENSMBL to insure that the ensemble is being constructed 
properly. 

2. Output Description 
Program SMUMEAN computes the non-stationary mean and 
writes an output file consisting of two columns. The first 


column is ensemble sample number (or row number, running from 


1 to LASTI). Knowing the sample rate, the time scale may be 
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inferred from sample number. The second column contains the 
non-stationary mean values. The following pages present the 


source code for program ENSEMBL. 
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PROGRAM SMUMEAN 
Fe TT I TTI RRR REE KK RKREREEKERRKEEEKERKKERE KEKE 


Program to perform low-pass FFT digital filter smoothing and then 
ensemble averaging on the smoothed turbulent boundary layer velocity 
data. In other words, this program computes the time varying mean of 
a non-stationary time series. 


The cutoff frequency of the low pass filter is specified by the 
size of the smoothing window, PTS, the sample rate and the number 
of points to be smoothed, M. See the table below: 


Window Size -3 dB Cutoff Frequency (FS=30000,M=2048) 


2 8120 
5 3250 
10 1625 
15 1080 
20 810 
25 650 
30 540 
a5 465 
40 405 
50 325 


For other sample rates and window sizes, the cutoff frequency may 
be determined from the equation: 


FC = (0.29289 * (FS/PTS) **2)**0.5 


where FC is the cutoff frequency, and FS is the sample rate. 

Note, the input array to the smoothing subroutine does not have to be 
an integer power of two because the subroutine zero pads as necessary 
to the next highest integer power of two greater than the length of 
the data vector (in this case the length of the columns of ENSMBL, 
which is LASTI.* 


The input data file is assumed to be in output format of the A/D 
conversion system, ie. a sequential file with alternating 

values of velocity (U) and propellor pulse (PROP), starting with 
velocity. Smoothing is accomplished with subroutine SMOOFT which 
performs frequency domain digital (low pass) filtering of 

the velocity signal. Program calculates ensemble mean velocity 
from the smoothed velocity values. Variance and turbulence 
intensity are computed in program TURBIN. Array size regqirements 
prevent doing both in one pass through the data. 


Key variables are defined below: 


ENSMBL = The primary data array which contains up to 50 pulses (columns) 
each of length up to 2048 points. The pulses are synchronized in 
time (an ensemble) such that each point ina pulse "lines up" 
with the others. 

U = The measured velocity as read from the input file into array ENSMBL. 

PROP = The propeller pulse timing marker signal. Also read from the input. 

IPLS = Number of pulses to be read from the input file. 

LASTI = Number of points in the columns of ENSMBL. (ie. rows in ENSMBL) 
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NPLS = (J-1) = Total number of pulses (columns in ENSMBL). 
AVG = The ensemble averaged then FFT smoothed value of U. (Output) 


Programmed by: Donald K. Johnson 
Sept. 1988,NPS 
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DIMENSION ENSMBL(2048,50) ,AVG(2048) 
DIMENSION UCOL(2048) 

INTEGER I, IPLS,LASTI,J,K,L 
CHARACTER*12 INFILE,OUTFILE 


c 
PRINT '(/A/)',' Low-Pass Filter Smoothing then Ensemble Average' 
Cc 
PRINT '(/A\)',' ENTER THE INPUT DATA FILE NAME: ! 
READ (*,'(A)') INFILE 
Cc 
Cc PRINT '(/A\)','* ENTER THE NUMBER OF PROP PULSES YOU WANT: '! 
e READ (*,*) IPLS 
IPLS=50 
C 
PRINT '{/A\)*,' ENTER THE OUTPUT DATA FILE NAME 
READ (*,'(A)') OUTFILE 
c 
PRINT ‘(/A\)*,' ENTER THE FFT SMOOTHING WINDOW SSIZE. 
READ(*,*) PTS 
S 
OPEN (10, FILE=INFILE,STATUS='OLD‘ , FORM='FORMATTED' ) 
OPEN (11, FILE=OUTFILE, STATUS='NEW') 
Cc 
LASTI=2048 
J=1 
I=0 
PRINT'(/A/)',' Smoothing then ensemble averaging:' 
PRINT'(/A/)',* Reading input file... 
cS 
10 CONTINUE 
wae Read the data into array "ENSMBL", each column is a pulse **#* 
READ (10,*,END=20,ERR=20) U, PROP 
L=1+1 
eet But skip the first partial pulse **#*#* 
IF (J .GT. 1) ENSMBL(1T,J-1)-U 
elalial If you have a new pulse, start a new column 
* prop trigger currently set to 1. volt *** 
IF(I «GE. 20 .AND. PROP .GT. 1. )@ THEN 
IF (J .EQ. 1) THEN 
WRITE(*,*)'Skipped first',I,' points (not a full pulse).' 
ELSE 
WRITE(*,*) ‘No. of points in column’ 73-1) 1s 7a 
ENDIF 
kkk Find the shortest column *** 


IF({I .LT. LASTI .AND. J .GT. 1) LAstien 
IF (J=-1 .GE. IPLS) Go To zs 
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J=J+1 
I=0 
ENDIF 
GO TO 10 
20 CONTINUE 
PRINT '(/A/)',' End of File or read error encountered - deletin 
lg last partial colunn....' 
J=J=-1 
28 CONTINUE 
XN=FLOAT (J-1) 


NPLS=J-1 
Cc 
WRITE(*,*) ' The final number of pulses is',J-1,' each of len 
1gth', LASTI 
Cc 


xxx Perform the FFT smoothing on each column of the array ENSMBL *** 
DO 300 L=1,NPLS 
WRITE(*,*) 'Smoothing column',L 
DO 310 I=1,LASTI 
UCOL(I) =ENSMBL(I,L) 


310 CONTINUE 
CALL SMOOFT(UCOL, LASTI, PTS) 
kee Write the smoothed velocity back into ENSMBL *** 


DO 320 I=1,LASTI 
ENSMBL(1I,L) =UCOL(T) 


320 CONTINUE 
300 CONTINUE 
c 
**x* Compute the ensemble mean *** 
c 
DO 22 I=1,LASTI 
SUM=0. 
SUMSQ=0. 


DO 30 L=1,NPLS 
SUM=SUM+ENSMBL(I,L) 
30 CONTINUE 
AVG (1) =SUM/XN 
22 CONTINUE 
xz*x Write the output ****** 
DO 40 I=1,LASTI 
XI=FLOAT (I) 
WRITE(11,200) XI,AVG(TI) 


200 FORMAT (‘F5.0,F8.1) 
40 CONTINUE 

CLOSE(10) 

CLOSE (11) 

END 


CR eR RRR RRR RRR RRR RRR RRR RRA AAREARR AEE 
SUBROUTINE SMOOFT(Y,N,PTS) 


Y IS THE DATA ARRAY, N IS THE LENGTH OF THE DATA ARRAY, AND PTS IS 
THE SMOOTHING WINDOW WIDTH. 


This subroutine is from the book “Numerical Reciepes" by Press, 
Flannery, Teukolsky and Vetterling, Pub. 1986, Cambridge Univ. Press 
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c 
C 


ANN 


ANAANNANANNA 


PARAMETER (MMAX=2048) 
DIMENSION Y(MMAX) 


M IS CALCULATED TO BE THE SMALLEST POWER OF 2 WHICH IS GREATER THAN 
THE NUMBER OF POINTS N. THE REMAINDER OF THE ARRAY IS ZERO PADDED. 


M=2 
NMIN=N+2. *PTS 
IF(M.LT.NMIN) THEN 
M=2*M 
GO TO 1 
ENDIF 
IF(M.GT.MMAX) PAUSE 'MMAX too small' 
CONST=(PTS/M) **2 
Y1=Y(1) 
YN=Y(N) 
RN1=1./(N-1.) 


REMOVE THE LINEAR TREND 


DO 11 J=1,N 
Y(J)=¥ (J) -RN1*(Y1* (N-J)+YN*(J-1)) 
1 CONTINUE 
IF (N+1.LE.M) THEN 


ZERO PAD 


DO 12 J=N+1,M 
Y(J)=0. 
2 CONTINUE 
ENDIF 


TRANSFORM TO FREQUENCY DOMAIN 


MO2=M/2 
CALL REALFT(Y,MO2,1) 
¥(1)=Y¥(1)/MO2 

FAC=1. 


WRITE(11,*) ' BIN (I) Y(F)! 
DO 99 I=1,M 
WRITE(11,*) I,Y(I) 
99 CONTINUE 


WINDOW FUNCTION ~ (USES THE PARABOLIC WELCH WINDOW) 


WRITE(11,*) ' WINDOW FUNCTION: ' 
DO 13 J=1,MO02-1 
K=2*J4+1 
IF(FAC.NE.0.) THEN 
FAC=AMAX1(0., (1--CONST*J**2) /MO2) 
WRITE(11,*) K,K+1,FAC 
Y (K) =FAC*Y (K) 
Y (K+1) =FAC*Y (K+1) 
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Cc 
c 
Cc 
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ELSE 
Y(K) =0. 
Y (K+1)=0. 
ENDIF 
CONTINUE 


HANDLES THE LAST POINT 


FAC=AMAX1(0., (1.-0.25*PTS**2) /MO2) 
Y (2) =FAC*Y(2) 

WRITE(11,*) ' WINDOW x DATA:' 

DO 98 I=1,M 

WRITE(11,*) I,Y¥(I) 


98 CONTINUE 


TRANSFORM BACK TO TIME DOMAIN 
GALL REALFI(Y,MO2,=21) 
RESTORE LINEAR TREND 


DO 14 J=1,N 
Y (J) =RN1* (Y1*(N-J) +YN*(J-1))+Y(J) 
CONTINUE 
RETURN 
END 
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SUBROUTINE REALFT(DATA,N,ISIGN) 


CALCULATES THE FOURIER TRANSFORM OF A SET OF 2*N REAL VALUED DATA 
POINTS (ARRAY 'DATA'). ISIGN IS 1 FOR THE FOURIER TRANSFORM (TIME 
TO FREQ. DOMAIN) AND ISIGN IS -1 FOR INVERSE TRANSFORMATIOM. 

THIS SUBROUTINE SETS UP THE DATA TO BE EFFICIENTLY TRANSFORMED BY 
SUBROUTINE ‘FOURI'. 


This subroutine is from the book "Numerical Reciepes" by Press, 
Flannery, Teukolsky and Vetterling, Pub. 1986, Cambridge Univ. Press 


REAL*8 WR,WI,WPR,WPI,WTEMP, THETA 
DIMENSION DATA(*) 

THETA=6 .. 28318530717959D0/2.0D0/DBLE(N) 
C1=0.5 


FORWARD TRANSFORM BY FOUR] 


TE (ISIGN.EQ.1) THEN 


C2=-0.5 
CALL FOUR1(DATA,N,+1) 
ELSE 


SET UP FOR INVERSE TRANSFORM 


C2=0.5 
THETA=-THETA 
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ENDIF 

WPR=-2.0D0*DSIN(0.5D0*THETA) **2 

WPI=DSIN (THETA) 

WR=1.0D0+WPR 

WI=WPI 

N2P3=2*N+3 

DO 11 I=2,N/2+1 
I1=2*I-1 
I2=I1+1 
I3=N2P3~-I2 
14=I3+1 
WRS=SNGL(WR) 
WIS=SNGL(WI) 
H1R=C1* (DATA(I1)+DATA(I3)) 
H1I=C1* (DATA(I2) -DATA(I4) ) 
H2R=-C2* (DATA (I2)+DATA(I4) ) 
H2I=C2* (DATA (I1) -DATA(I3) ) 
DATA (11) =H1R+WRS*H2R-WIS*H2I 
DATA (I2) =H1I+WRS*H2I+WIS*H2R 
DATA (13) =H1R-WRS *H2R+WIS*H2I 
DATA (I4)=-H1I+WRS *H2I+WIS*H2R 
WIEMP=WR 
WR=WR*WPR-WI *WPI+WR 
WI=WI*WPR+WTEMP*WPI+WI 

sal CONTINUE 

IF (ISIGN.EQ.1) THEN 
H1R=DATA(1) 
DATA (1) =H1R+DATA (2) 
DATA (2) =H1R-DATA(2) 

ELSE 
H1R=DATA(1) 
DATA(1)=C1* (H1R+DATA (2) ) 
DATA (2) =C1* (H1R-DATA (2) ) 


CG 
C INVERSE TRANSFORM 
c 
CALL FOUR1(DATA,N,-1) 
ENDIF 
RETURN 
END 
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SUBROUTINE FOUR1(DATA,NN,ISIGN) 
G wee e@eee eee eee eo wwe eee eee eee eee eee ww OO ee eee ewe we ee we we eM ew ew we wee ew eee 
C COMPUTES THE COMPLEX DFT OR INVERSE DFT OF ARRAY 'DATA‘' OF 
C LENGTH 'NN‘'. ISIGN IS 1 FOR THE DFT AND -1 FOR THE INVERSE 
CG. DEI. 
e 
C This subroutine is from the book "Numerical Reciepes" by Press, 
C Flannery, Teukolsky and Vetterling, Pub. 1986, Cambridge Univ. 
Ee See wee ewe e@ e@ Oe Be SS ew @ @ @ @ & & or ww oe ee © Oe SO ee 8 Oe 0 

REAL*8 WR,WI,WPR,WPI,WTEMP, THETA 

DIMENSION DATA(*) 

N=2*NN 

J=1 

DO 11 I=1,N, 2 
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[P(J.Gl. 1) THEN 
TEMPR=DATA (J) 
TEMPI=DATA(J+1) 
DATA (J) =DATA (TI) 
DATA (J+1) =DATA(I+1) 
DATA (I)=TEMPR 
DATA (I+1)=TEMPI 
ENDIF 
M=N/2 
IF ((M.GE.2).AND.(J.GT.M)) THEN 
J=J-M 
M=M/2 
GO TO l 
ENDIF 
J=J+M 
CONTINUE 
MMAX=2 
IF (N.GT.MMAX) THEN 
ISTEP=2 *MMAX 
THETA=6 . 28318530717959D0/ (ISIGN*MMAX) 
WPR=-2.D0*DSIN(0O.5DO*THETA) **2 
WPI=DSIN (THETA) 
WR=1.D0 
WI=0.D0 
DO 13 M=1,MMAX,2 
DO 12 I=M,N,ISTEP 
J=I+MMAX 
TEMPR=SNGL (WR) *DATA(J) -SNGL(WI) *DATA(J+1) 
TEMPI=SNGL(WR) *DATA(J+1)+SNGL(WI) *DATA(J) 
DATA (J) =DATA(I)-TEMPR 
DATA (J+1) =DATA(I+1) -TEMPI 
DATA(I)=DATA(I)+TEMPR 
DATA (I+1) =DATA(I+1)+TEMPI 
CONTINUE 
WTEMP=WR 
WR=WR*WPR-WI*WPI+WR 
WI=WI*WPR+WTEMP*WPI+WI 
CONTINUE 
MMAX=ISTEP 
GO TO 2 
ENDIF 
RETURN 
END 
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E. USER’S GUIDE FOR PROGRAM MEANSMU 


1. Input Instructions 


¢ Insure that the executable file for program MEANSMU 
(MEANSMU.EXE) and the input data file are in the same 
subdirectory of a disk drive with sufficient space 
remaining to write the output file. 


¢* Initiate program MEANSMU (type "MEANSMU" at the DOS 
prompt. 


¢- The program will prompt the user for the name of the 
input data file. Enter the full name of the input data 
file including the extension. 


- The program will prompt the user for the name of the 
output file which will be created by the program. Enter 
the desired file name. Do not use the name of an 
existing file or the program will terminate. 


- The program will prompt the user for the FFT smoothing 
window size. The window size specifies the cutoff 
frequency of the low pass frequency domain filter. The 
low pass filtered data produces an interim "local" mean 
or local smoothed velocity which is then ensemble 
averaged. Table IV-2 in Chapter IV presents some 
examples of the relationship between moving average 
window size and its associated low pass filter cutoff 
frequency. Other cutoff frequencies may be calculated 
uSing equation A-11 in Appendix A. 


- The program will then read the input data file, compute 
the non-stationary mean and write the output data file. 
Remember to observe the column lengths of the array 
ENSMBL to insure that the ensemble is being constructed 
properly. 

2. Output Description 
Program MEANSMU computes the non-stationary mean and 
writes an output file consisting of two columns. The first 


column is ensemble sample number (or row number, running from 


1 to LASTI). Knowing the sample rate, the time scale may be 
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inferred from sample number. The second column contains the 
non-stationary mean values. The following pages present the 


source code for program ENSEMBL. 
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PROGRAM MEANSMU 
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Program to perform ensemble averaging then frequency domain low-pass 
digital filter smoothing on an ensemble of propeller pulses. 

The data is ensemble averaged first, then the average is filtered. 
In other words, this program computes the time varying mean of a 
non-stationary time series. 


The cutoff frequency of the low pass filter is specified by the 
Size of the smoothing window, PTS, the sample rate and the number 
of points to be smoothed, M. See the table below: 


Window Size -3 dB Cutoff Frequency (FS=30000,M=2048) 
2 8120 
5 3250 

10 1625 
iS 1080 
20 810 
Zo 650 
30 540 
is. 465 
40 405 
50 323 


For other sample rates and window sizes, the cutoff frequency may 
be determined from the equation: 


FE = (0.29289 * (ES/PITS) **2) "2075 


where FC is the cutoff frequency, and FS is the sample rate. 

Note, the input array to the smoothing subroutine does not have to be 
an integer power of two because the subroutine zero pads as necessary 
to the next highest integer power of two greater than the length of 
the data vector (in this case the length of the columns of ENSMBL, 


which is LASTI. 


Input data file is assumed to be in output format of the A/D 
conversion system, ie. a sequential ASCII file with alternating 
values of velocity (U) and propeller pulse (PROP), starting with 
velocity. Smoothing is accomplished with subroutine SMOOFT 

which performs frequency domain digital (low pass) filtering of 

the velocity signal. Program also calculates variance for the 
ensemble using the ensemble averaged-then-smoothed values of velocity. 


Key variables are defined below: 


ENSMBL = The primary data array which contains up to 50 pulses (columns) 
each of length up to 2048 points. The pulses are synchronized in 
time (an ensemble) such that each point ina pulse "lines up" 
with the others. 

U = The measured velocity as read from the input file into array ENSMBL. 

PROP = The propeller pulse timing marker signal. Also read from the input. 

IPLS = Number of pulses to be read from the input file. 

LASTI = Number of points in the columns of ENSMBL. (ie. rows in ENSMBL) 

NPLS = (J-1) = Total number of pulses (columns in ENSMBL). 
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AVG = The ensemble averaged then FFT smoothed value of U. (Output) 

VAR = Ensemble statistical variance of U. (Output) 

TURBIN = "Pseudo" ensemble turbulence intensity of U. (Not based on the 
local smoothed velocity but the final mean velocity. Similar to 
the statistical variance.) Output 

UINF = Free stream velocity external to the boundary layer 

PTS = Size of the smoothing window 


Programmed by: Donald K. Johnson 
July 1988,NPS 


TO OECECEEPCECE CESSES ESLER ERE SESE LER REE ESLER ERE LEE ESE LES SEES CCR SALES 2 oe 


DIMENSION ENSMBL(2048,50) ,AVG(2048) , VAR(2048) , TURBIN (2048) 
INTEGER I,IPLS,LASTI,J,K,L 

REAL UINF 

CHARACTER*12 INFILE, OUTFILE 


PRINT '(/A/)',' Ensemble Averaging and Low-Pass Filter Smoothing' 
PRINT '(/A\)',' ENTER THE INPUT DATA FILE NAME: ' 

READ (*,'({A)') INFILE 

PRINT '(/A\)',' ENTER THE NUMBER OF PROP PULSES YOU WANT: ' 

READ (*,*) IPLS 

IPLS=50 

PRinNt ) {(7A\)',° ENTER THE OUTPUT DATA FILE NAME: ‘ 


READ (*,'(A)') OUTFILE 


pseudo turbulence intensity calcs currently inactive *** 

PRINT '(/A)',' ENTER U-INFINITY (REF. VELOCITY FOR TURB.INTENSITY 
1 CALCS) - a REAL number: ! 

READ(*,*) UINF 

UINF=UINF*1. 

UINF=1. 


PRINT'(/A\)',' ENTER THE FFT SMOOTHING WINDOW SIZE: ' 
READ(*,*) PTS 


OPEN (10, FILE=INFILE, STATUS='OLD' , FORM=' FORMATTED ' ) 
OPEN(11,FILE=OUTFILE, STATUS='NEW' ) 


LASTI=2048 

J=1 

I=0 

PRINT'(/A/)',' Computing ensemble then FFT smoothed average:' 
PRIND (7/7) ,.' READING INPUT FILE....' 


CONTINUE 
Read the data into array “ENSMBL", each column is a pulse *** 
READ (10, *,END=20,ERR=20) U,PROP 
I=I+1 
But skip the first partial pulse ***x*** 
IF (J .GT. 1) ENSMBL(I,J-1) =U 
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eae If you have a new pulse, start a new column (prop trigger currently 
* set "£0 1. Volt) 
IF(I .GE. 20 .AND. PROP -GT. 2.) THEN 
IF (J .EQ. I), THEN 
WRITE(*,*) 'Skipped first',I,' points (not a full pulse).' 


ELSE 
WRITE(*,*) 'No. of points in column',J-1,' is',I 
ENDIF 
kat Find the shortest column *** 


IF(I .LT. LASTI .AND. J .GT. 1) LASTI=I 
IF (J-1 .GE. IPLS) GO TO 28 
J=J+1 
I=0 
ENDIF 
GO TO 10 
20 CONTINUE 
WRITE(*,'(/A/)')' End of File or read error encountered - delet 
ling last partial column... 
J=J-1 
28 CONTINUE 
XN=FLOAT (J-1) 
NPLS=J-1 


WRITE(*,*)'The final number of pulses is',J-1,' each of length', 

1 LASTI 
***x Compute the ensemble mean, variance and “pseudo” turbulence intensity 
* based on user input value of U - infinity (free stream velocity) 
* and the smoothed mean velocity. 


*x*x compute ensemble mean *** 
DO 22 I=1,LASTI 
SUM=0O. 
DO 30 L=1,NPLS 
SUM=SUM+ENSMBL(I,L) 
30 CONTINUE 
AVG (I) =SUM/XN 
22 CONTINUE 
*x* smooth the ensemble average *** 
PRINT'(/A/)',' Smoothing ensemble average....' 
CALL SMOOFT (AVG, LASTI, PTS) 
z*kx compute variance and pseudo turbulence inten. *** 
PRINT'(/A/)',' Computing turbulence intensity....' 
DO 50 I=1,LASTI 
SUMSQ=0. 
DO 52 L=1,NPLS 
SUMSQ=SUMSQ+ (ENSMBL(I, L) -AVG(I)) **2. 
52 CONTINUE 
VAR(I) = SUMSQ/(XN-1.) 
TURBIN(I) = ((SUMSQ/XN) **.5) /UINF 
50 CONTINUE 
z*z*x Write the output ****** 
PRINT "(74/7 ))).,' OUGpucting 
DO 40 I=1,LASTI 
XI=FLOAT (TI) 
eat note: variance and pseudo turbulence intensity is currently 
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not outpucs*** 
WRITE(11,200) XI,AVG(I) 


200 FORMAT (F5.0,F8.1) 
40 CONTINUE 

CLOSE(10) 

CLOSE (11) 

END 


CRRRAAEKAAAEAEAKAAEKAEKKERKERAEEKREKEKEEKREKEAERERAEKAKEREKKEKREKEKRRKEKRKEKKRKEKREKKKAK 
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SUBROUTINE SMOOFT(Y,N, PTS) 


Y IS THE DATA ARRAY, N IS THE LENGTH OF THE DATA ARRAY, AND PTS IS 
THE SMOOTHING WINDOW WIDTH 


This subroutine is from the book "Numerical Reciepes" by Press, 
Flannery, Teukolsky and Vetterling, Pub. 1986, Cambridge Univ. Press 


PARAMETER (MMAX=2048) 
DIMENSION Y(MMAX) 


M IS CALCULATED TO BE THE SMALLEST POWER OF 2 WHICH IS GREATER THAN 
THE NUMBER OF POINTS N. THE REMAINDER OF THE ARRAY IS ZERO PADDED. 


M=2 
NMIN=N+2.*PTS 
IF(M.LT.NMIN) THEN 
M=2*M 
elo) bee al 
ENDIF 
IF(M.GT.MMAX) PAUSE 'MMAX too small' 
CONST=(PTS/M) **2 
Y1=¥(1) 
YN=Y (N) 
RN1=1./(N-1.) 


REMOVE THE LINEAR TREND 


Doe li csi N 
¥ (J) =Y (J) -RN1* (Y1*(N-J)+YN*(J-1) ) 

CONTINUE 

IF (N+1.LE.M) THEN 


ZERO PAD 


DO 12 J=N+1,M 
Y(J)=0. 
CONTINUE 
ENDIF 


C TRANSFORM TO FREQUENCY DOMAIN 


Cc 


MO2=M/2 

CALL REALFT(Y,MO2,1) 
¥(1)=¥ (1) /MO2 

FAC=1. 


eS 


WRITE(11,*)' BIN (I) Y(F)! 
DO 99 I=1,M 
WRITE(11,*) I,Y(I) 
99 CONTINUE 


WINDOW FUNCTION - (USES THE PARABOLIC WELCH WINDOW) 


ANANAAANNA 


WRITE(11,*) ' WINDOW FUNCTION:' 
DO 13 J=1,M02-1 
K=2*J+1 
IF(FAC.NE.0.) THEN 
FAC=AMAX1(0., (1.-CONST*J**2) /MO2) 
Cc WRITE(11,*) K,K+1,FAC 
Y (K) =FAC*Y (K) 
Y (K+1) =FAC*Y (K+1) 
ELSE 
Y (K) =0. 
Y (K+1) =0. 
ENDIF 
3 CONTINUE 


HANDLES THE LAST POINT 


ANNRr 


FAC=AMAX1(0., (1.-0.25*PTS**2) /MO2) 
Y (2) =FAC*Y(2) 
WRITE(11,*) ' WINDOW x DATA:' 
DO 98 I=1,M 
WRITE(11,*) I,Y(I) 
98 CONTINUE 


TRANSFORM BACK TO TIME DOMAIN 


ANAAANNAN 


CALL REALFT(Y,MO2,-1) 


RESTORE LINEAR TREND 


aan 


DO 14 J=1,N 
Y (J) =RN1* (Y1* (N-J)+YN*(J-1) )+Y(J) 
14 CONTINUE 
RETURN 
END 
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SUBROUTINE REALFT(DATA,N,ISIGN) 
CALCULATES THE FOURIER TRANSFORM OF A SET OF 2*N REAL VALUED DATA 
POINTS (ARRAY 'DATA'). ISIGN IS 1 FOR THE FOURIER TRANSFORM (TIME 
TO FREQ. DOMAIN) AND ISIGN IS -1 FOR INVERSE TRANSFORMATIOM. 
THIS SUBROUTINE SETS UP THE DATA TO BE EFFICIENTLY TRANSFORMED BY 
SUBROUTINE 'FOURI1'. 


This subroutine is from the book "Numerical Reciepes" by Press, 
Flannery, Teukolsky and Vetterling, Pub. 1986, Cambridge Univ. Press 
REAL*8 WR,WI,WPR,WPI,WTEMP, THETA 
DIMENSION DATA(*) 


AANANANANNANN 
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THETA=6. 28318530717959D0/2.0D0/DBLE(N) 
C1=0.5 


FORWARD TRANSFORM BY FOUR] 


aaa 


IF (ISIGN.EQ.1) THEN 
C2=-0.5 
CALL FOUR1(DATA,N,+1) 
ELSE 


SET UP FOR INVERSE TRANSFORM 


Qa0aAN 


C2=0.5 
THETA=-THETA 

ENDIF 

WPR=-2.0D0*DSIN(0.5DO*THETA) **2 

WPI=DSIN (THETA) 

WR=1.0D0+WPR 

WI=WPI 

N2P3=2*N+3 

DO 11 I=2,N/2+1 
I1=2*I-1 
I2=I1+1 
I3=N2P3-I2 
I4=I3+1 
WRS=SNGL (WR) 
WIS=SNGL(WI) 
H1R=C1* (DATA(I1)+DATA(I3)) 
H1I=C1* (DATA(I2) -DATA(I4) ) 
H2R=-C2* (DATA(I2)+DATA(I4) ) 
H2I=C2* (DATA(I1)-DATA(I3) ) 
DATA (I1) =HLR+WRS*H2R-WIS*H2I 
DATA (12) =H1I+WRS*H21+WIS*H2R 
DATA (13) =H1R-WRS*H2R+WIS*H2I 
DATA (14) =-H1I+WRS*H2I+WIS*H2R 
WTEMP=WR 
WR=WR*WPR-WI1*WPI+WR 
WI=WI*WPR+WTEMP*WPI+WI 

ital CONTINUE 

IF (ISIGN.EQ.1) THEN 
H1R=DATA (1) 
DATA (1) =H1R+DATA(2) 
DATA (2) =H1R-DATA (2) 

ELSE 
H1R=DATA (1) 
DATA (1) =C1* (H1R+DATA (2) ) 
DATA (2) =C1* (H1R-DATA (2) ) 


e 
C INVERSE TRANSFORM 
e 
CALL FOUR1(DATA,N,-1) 
ENDIF 
RETURN 
END 
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SUBROUTINE FOUR1(DATA,NN, ISIGN) 


Crm mw wm mm mm mmm ew ww we wwe ew ww www www wee re tr ene ew eer err errr rere 
C COMPUTES THE COMPLEX DFT OR INVERSE DFT OF ARRAY 'DATA' OF 
C LENGTH 'NN'. ISIGN IS 1 FOR THE DFT AND =-1 FOR THE INVERSE 
CG eDET. 
c 
C This subroutine is from the book "Numerical Reciepes" by Press, 
C Flannery, Teukolsky and Vetterling, Pub. 1986, Cambridge Univ. Press 
Cw ree em mn wo an www ww ww www ew we ee ew ew www wm mm em me mm em em mee eee err ore rerrerrrrr-=- 
REAL*8 WR,WI,WPR,WPI,WTEMP, THETA 
DIMENSION DATA(*) 
N=2*NN 
J=1 
DO 11 I=1,N,2 
LF (9.61.1) THEN 
TEMPR=DATA (J) 
TEMPI=DATA(J+1) 
DATA (J) =DATA (TI) 
DATA (J+1)=DATA(I+1) 
DATA(I)=TEMPR 
DATA (I+1) =TEMPI 
ENDEF 
M=N/2 
nt IF ((M.GE.2).AND.(J.GT.M)) THEN 
J=J-M 
M=M/2 
GO TO l 
ENDIF 
J=J+M 
peat CONTINUE 
MMAX=2 
2 IF (N.GT.MMAX) THEN 
ISTEP=2 *MMAX 
THETA=6.28318530717959D0/ (ISIGN *MMAX) 
WPR=-2.D0*DSIN(0.5DO*THETA) **2 
WPI=DSIN (THETA) 
WR=1.D0 
WI=0.D0 
DO 13 M=1,MMAX,2 
DO 12 1-M iN, LSteEP 
J=I+MMAX 
TEMPR=SNGL(WR) *DATA(J) -SNGL(WI) *DATA(J+1) 
TEMPI=SNGL(WR) *DATA(J+1)+SNGL(WI) *DATA(J) 
DATA(J)=DATA(I)-TEMPR 
DATA(J+1)=DATA(I+1) -TEMPI 
DATA(I) =DATA(I)+TEMPR 
DATA (I+1) =DATA(I+1)+TEMPI 
Le CONTINUE 
WTEMP=WR 
WR=WR*WPR-WI*WPI+WR 
WI=WI*WPR+WTEMP*WPI+WI 
13 CONTINUE 
MMAX=ISTEP 
GO TO 2 
ENDIF 
RETURN 
END 


GC 


F. 


USER’S GUIDE FOR PROGRAM TURINTMA 


1. Input Instructions 


Insure that the executable file for program TURINTMA 
(TURINTMA.EXE) and the input data file are in the same 
subdirectory of a disk drive with sufficient space 
remaining to write the output file. 


Initiate program TURINTMA (type "TURINTMA" at the DOS 
prompt). 


The program will prompt the user for the name of the 
input data file. Enter the full name of the input data 
file including the extension. 


The program will prompt the user for the name of the 
output file which will be created by the program. Enter 
the desired file name. Do not use the name of an 
existing file or the program will terminate. 


The program will prompt the user for the moving average 
window size. The window size specifies the number of 
samples which are time averaged to produce an interim 
"local" mean or local smoothed velocity which is used to 
compute the fluctuating velocity component. The window 
Size should be an odd number. Table IV-12 in Chapter IV 
presents some examples of the relationship between moving 
average window size and its associated low pass filter 
cutoff frequency. Other cutoff frequencies may be 
determined by solving equation A-7 in Appendix A. 
Program FCMA is provided to solve equation A-7. 


The program will prompt the user for the reference 
velocity used to normalize the turbulence intensity 
values. The mean velocity at the outer edge of the 
boundary layer (U,,,,) 1s the typical reference velocity. 
The user will need to determine this value PELOT > CO 
executing this program. 


The program will then read the input data file, compute 
the turbulence intensity and write the output data file. 
Remember to observe the column lengths of the array 
ENSMBL to insure that the ensemble is being constructed 
properly. 


1 (31 7 


2. Output Description 

Program TURINTMA computes the non-stationary 
turbulence intensity, then writes an output file consisting 
of two columns. The first column is ensemble sample number 
(or row number, running from 1 to LASTI). Knowing the sample 
rate, the time scale may be inferred from sample number. The 
second column contains the non-stationary turbulence intensity 
values. The following pages present the source code for 


program ENSEMBL. 
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$DEBUG 


PROGRAM TINTMA 


CRRA RARRAKRERRREEKEREKEKEKERAERKEEKEEEKEERERAEEAERAAEREKEKEKEREKAREKKERER 


Program computes turbulence intensity as the square root of the ensemble 
average of the mean square fluctuating velocity component. The 

fluctuating velocity component is the difference between the 

instantaneous velocity and the local smoothed velocity, computed 

using a moving average smoothing algorithm. The turbulence intensity is 
normalized using an input reference velocity. Note that the moving average 
operation sets the first (L-1)/2 points and the last (L=-1)/2 points 

of the mean velocity output to zero (where L is the moving average 

window size). 


The relationship between smoothing window size and the corresponding 
cutoff frequency of the smoothing operation is given by: 


SIN(( (uel) /2 +e. > ees PI*FC/FS) / L*SIN(PI®FC/FS) = 027071 


where FC is the cutoff frequency, FS is the sample rate and L is the 
odd numbered window size. 


Input data file is assumed to be in output format of the A/D 
system, i.e. a sequential file with alternating values 

of velocity (U) and propeller pulse (PROP), starting with velocity. 
The program first separates the single time history file into an 
ensemble of pulses using the value of the propeller trigger (PROP) 
as the marker for the beginning of a pulse. 


Key variables are defined below: 


ENSMBL = The primary data array which contains up to 50 pulses (columns) 
each of length up to 2048 points. The pulses are synchronized in 
time (an ensemble) such that the beginning of each column (i.e. 
the first row) represents the beginning of each pulse. 

U = The measured velocity as read from the input file into array ENSMBL. 

PROP The propeller pulse timing marker signal. Also read from the input. 

IPLS Number of pulses to be read from the input file. 

LASTI = Number of points in the columns of ENSMBL. (ie. rows in ENSMBL) 

NPLS = (J-1) = Final total number of pulses (columns in ENSMBL). 

L = The width of the moving average window (should be an odd number). 

A(I) = The current column of ENSMBL to be smoothed using subroutine MOVAVE 

B(I) = The current smoothed column of ENSMBL. (Output of MOVAVE) 

TURBIN(I) =:The average turbulence intensity for NPLS propeller pulses. 

UEDGE = The reference velocity used to normalize the turbulence intensity 
(the mean velocity at the outer edge of the boundary layer is 
usually used - input from keyboard). 


Programmed by: Donald K. Johnson 
AUGUST, 1988,NPS 
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PARAMETER (MMAX=2048) 
DIMENSION ENSMBL(2048,50) ,A(MMAX) , B(MMAX) , TURBIN (MMAX) 
INTEGER L,NPTS,ILAST,NPLS,IPLS 

REAL UEDGE 
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CHARACTER#1 2 INFI LE ,OUTFIEE 
PRINT'(/A/)',' Turbulence Intensity from Moving Average Smoothing' 


PRINT '(/A\)',' ENTER THE INPUT DATA FILE NAME: '! 
READ (*,'(A)') INFILE 


PRINT '(/A\)',' ENTER THE OUTPUT FILE NAME: '! 
READ (*,'(A)') OUTFILE 


PRINT '(/A\)',' ENTER THE MOVING AVERAGE WINDOW SIZE" (AN ODD NUMB 
TER) 
READ (*,*) L 


PRINT '(/A\)',' ENTER U-EDGE (REFERENCE VELOCITY FOR TURBULENCE IN 
1TENSITY NORMALIZATION): '! 

READ(*,*) UEDGE 

UEDGE=UEDGE*1. 


OPEN (UNIT=10, FILE=INFILE, STATUS='OLD') 
OPEN (UNIT=11, FILE=OUTFILE, STATUS='NEW') 


ee separate the time history data into an ensemble of pulses *** 


LASTI=2048 
IPLS=50 
J=1 
I=0 
PRINT'(/A/)*,' READING INPUT FILE... -'! 
PRINT'(/A/)',' WARNING, Watch the number of points in each column 
1... they should all be about the same length.' 
x 
10 CONTINUE 
Lala! Read the data into array "ENSMBL", each column is a pulse *** 
READ (10,*,END=20,ERR=20) U,PROP 
I=I+1 
ake But skip the first partial pulse *****x* 
IF (J .GT. 1) ENSMBL(I,J-1)=U 
alata If you have a new pulse, start a new column. *** 
* Note, prop trigger is set to 1. volt to start a new pulse *** 


IF(I 3GE. 20 <ANDO.. PROP Gl.) oertneN 
TPC) EO.) THEN 
WRITE(*,*) 'Skipped first',I,' points (not a full pulse).' 


ELSE: 
WRITE(*,*)'No. of points in column’ 7J—17) 91-2 
ENDIF 
xe tt Find the shortest column *** 
IF(I .LT. LASTI .AND. J <GT. 1). EASsTti=n 
eae Quit loop if you have enough pulses *** 
IF (J=1 (GE. FPES)\/Go ro 28 
J=J+1 
I=0 
ENDIF 
GO TO 10 


20 CONTINUE 
WRITE(*,'(/A/)')' End of File encountered - deleting last parti 
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60 
eee 


50 
ket 


100 
70 


Crs 


* 


Rai CO umn... 
J=J-1 
CONTINUE 
XN=FLOAT (J-1) 
NPLS=J-1 


WRITE(*,*)'The final number of pulses is',NPLS,' each of length', 


1 LASTI 


PRINT '(/A/)',' Working...' 


apply the moving average window to each column (pulse) in array ENSMBL 


start the pulse loop *** 

DO 25 K=1,NPLS 
and start the time loop *** 
DO 15 I=1,LASTI 


copy the current (Kth) pulse into A for call to subroutine *** 


A(I)=ENSMBL(I, K) 
CONTINUE 
smooth the Kth pulse *** 
CALL MOVAVE(A,LASTI,L,B) 


compute the square of the fluctuating velocity component 
and write value back into ENSMBL 
DO 16 I=1,LASTI 
ENSMBL(1I,K)=(ENSMBL(I,K)-B(I)) **2. 
CONTINUE 
CONTINUE 


now ensemble average (across rows) to get the mean square fluctuating 


velocity. 
start with a time loop *** 
DO 50 I=1,LASTI 
now sum across each row ***% 
SUMU=0. 
DO 60 K=1,NPLS 
SUMU=SUMU+ENSMBL(I,K) 
CONTINUE 
compute the normalized turbulence intensity (the root mean square) 
TURBIN (I) =(SUMU/NPLS) **.5/UEDGE 
IF(I .LE. (L-1)/2 .OR. I .GT. LASTI-(L-1)/2 ) TURBIN(I)=0. 
CONTINUE 
write the output *** 
DO 70 I=1,LASTI 
XI=FLOAT (TI) 
WRITE(11,100) XI,TURBIN(I) 
FORMAT (1%,F5.0,1X%,F10. 5) 
CONTINUE 
CLOSE(10) 
SLOSE(11) 
END 
RaEAKAERAAEKAKAAKKKEKKKKKEKKEKKKRKKEKKKKKaAKKEKKKEKKKKKKKEKKKKEKKEKRKRKKKEKKKaAKKEKRE 


SUBROUTINE MOVAVE(A,N,L,B) 
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Program performs a simple moving average of a time series to determine a 
local mean (smoothed) velocity. Note that the first L/2 values and the 
last L/2 values of the input time series are set to zero. Briefly, the 
moving average accomplishes a crude low pass filter type of smoothing 

on the input time series as follows. Starting with the first L (usually 
odd) points in the time series, a time average is computed for the L 
points. This average value becomes the smoothed velocity at the center 

of the L points, point (L-1)/2 for odd L. Then the moving average window 
is moved over one point and a new average is computed for the window. This 
is repeated until the Lth point of the window reaches the end of the time 
series. 


The key variables are defined below: 


A = A column array containing the time series to be averaged. 
N = The number of values in A. 
L = The size of the averaging window (an odd number is best because 


then the average of the small time segment is centered on the 
mid-point of the segment). Note, a larger window size corresponds to 
a lower low pass filter cutoff frequency. See the table below for the 
relation between window size and cutoff frequency. 

B = A column array containing the smoothed (averaged) time series values. 


DIMENSION A(1),B(1) 


AL=L 
MH=L/2 
zero out the answer array *** 
DO 18 I=1,N 
18 B(I)=0. 
K=L-1 
if L is .GT. N then quit *** 


EP CUNa ib) el,2,2 
2 NN=N-K 
start the main averaging loop, working on L/2 points on either 
side of the KK th data point. (KK is calculated below) 
DO 4 J=1,NN 
JN=J+K 
reset the accumulator *** 
S=0. 
sum across L points centered on point KK. 
DO 3 I=J,JN 
S=S+A(I) 
3 CONTINUE 
now divide by L to get the average of segment centered on point KK. 
KK=J +MH 
B (KK) =S/AL 
4 CONTINUE 
if L is even then take the average of two adjacent terms to get the value 
of the average velocity centered on the KK th point. 
IP etb=(27 2) 42) 1,6,24 
6 MS=MH+1 
MD=MH+NN-1 
DO 7 I=MS,MD 
TEMP=(B(I)+B(I+1))/2. 
B(I)=TEMP 


7 CONTINUE 
1 RETURN 
END 
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G. 


USER’S GUIDE FOR PROGRAM TURBIN 


1. Input Instructions 


Insure that the executable file for program TURBIN 
(TURBIN.EXE) and the input data file are in the same 
subdirectory of a disk drive with sufficient space 
remaining to write the output file. 


Initiate program TURBIN (type "TURBIN" at the DOS 
prompt). 


The program will prompt the user for the name of the 
input data file. Enter the full name of the input data 
file including the extension. 


The program will prompt the user for the name of the 
output file which will be created by the program. Enter 
the desired file name. Do not use the name of an 
existing file or the program will terminate. 


The program will prompt the user for the FFT smoothing 
window size. The window size specifies the cutoff 
frequency of the low pass frequency domain filter. The 
low pass filtered data produces an interim "local" mean 
or local smoothed velocity which is used to compute the 
fluctuating velocity component. Table IV-23 in Chapter 
IV presents some examples of the relationship between 
moving average window size and its associated low pass 
filter cutoff frequency. Other cutoff frequencies may 
be calculated using equation A-11 in Appendix A. 


The program will prompt the user for the reference 
velocity used to normalize the turbulence intensity 
values. The mean velocity at the outer edge of the 
boundary layer (U,,,,) is the typical reference velocity. 
The user will need to determine this value prior to 
executing this program. 


The program will then read the input data file, compute 
the turbulence intensity and write the output data file. 
Remember to observe the column lengths of the array 
ENSMBL to insure that the ensemble is being constructed 
properly. 
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2. Output Description 

Program TURBIN computes the non-stationary turbulence 
intensity, then writes an output file consisting of two 
columns. The first column is ensemble sample number (or row 
number, running from 1 to LASTI). Knowing the sample rate, 
the time scale may be inferred from sample number. The second 
column contains the non-stationary turbulence intensity 
values. The following pages present the source code for 


program ENSEMBL. 
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PROGRAM TURINT 
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Program computes the non-stationary turbulence intensity of an ensemble 
of propeller pulses by calculating the RMS of the fluctuating velocity. 
The mean in this case is an ensemble average. The fluctuating component 
is the difference between the instantaneous velocity and the local 
smoothed velocity, computed using a frequency domain low pass filter 
smoothing routine. The turbulence intensity is normalized using an 
input reference velocity. This program 1s the turbulence intensity 
counterpart to the non-stationary mean program - MEANSMU. 


The cutoff frequency of the low pass smoothing filter is specified by the 
size of the smoothing window, PTS, the sample rate and the number 
of points to be smoothed, M. See the table below: 


Window Size -3 dB Cutoff Frequency (FS=30000,M=2048) 
2 8120 
5 3250 

10 1625 
15 1080 
20 810 
25 650 
30 540 
35 465 
40 405 
50 S25 


For other sample rates and window sizes, the cutoff frequency may 
be determined from the equation: 


FC = (0.29289 * (FS/PTS) **2)**0.5 


where FC is the cutoff frequency, and FS is the sample rate. 

Note, the input array to the smoothing subroutine does not have to be 
an integer power of two because the subroutine zero pads as necessary 
to the next highest integer power of two greater than the length of 
the data vector (in this case the length of the columns of ENSMBL, 
which is LASTI. 


Input data file is assumed to be in output format of the Labtech 
Notebook system supporting the D/A board, ie. a sequential 

file with alternating values of velocity (U) and propellor pulse 
(PROP), starting with velocity. 


Key variables are defined below: 


ENSMBL = The primary data array which contains up to 50 pulses (columns) 
each of length up to 2048 points. The pulses are synchronized in 
time (an ensemble) such that each point ina pulse "lines up" 
with the corresponding points in other pulses. 

U = The measured velocity as read from the input file into array ENSMBL. 

PROP = The propeller pulse timing marker signal. Also read from the input. 

IPLS = Number of pulses to be read from the input file. 

LASTI = Number of points in the columns of ENSMBL. (ie. rows in ENSMBL) 

NPLS = (J-1) = Total number of pulses (columns in ENSMBL). 
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VAR = Ensemble variance of U. 
TURBIN = Ensemble turbulence intensity of U. 
UINF = Free stream velocity external to the boundary layer 


Programmed by: Donald K. Johnson 
July 1988,NPS 


RAARAARAEEAACRARKEKAAEKEKAAAKAKAEEKAKAEEEEEKREAKEAEEARKEEKAKEKRKAKRKKEKREKEKKEKEKKREKAKS 


DIMENSION ENSMBL(2048,50) ,VAR(2048) , TURBIN (2048) 
DIMENSION UCOL( 2048) 

INTEGER I, TPES; ASTI OFA, 

CHARACTER*#12 INFILE,OUTFILE 


PRINT ‘'(/A/)',‘ Turbulence Intensity from Frequency Domain Low Pas 
ls Filter Smoothing ' 
PRINT *(/A\)*,° ENTER THE INPUT DATA FILE NAME. 
READ (*,'(A)*) INFILE 
PRINT ‘(/A\)',' ENTER THE NUMBER OF PROP PULSES YOU WANT: ' 
READ (*,*) IPLS 
IPLS=50 
PRINT *(/A\)"," ENTER THE OUTPUT DATA PILE NAME 


Z 


ul 


1 


READ) (*, (A) *) ‘OUTEILE 


PRINT '(/A)‘,' ENTER U-EDGE (REFERENCE VELOCITY FOR TURBULENCE 
ENSITY NORMALIZATION): ' 

READ(*,*) UINF 

UINF=UINF#?1. 


PRINT'(/A\)',' ENTER THE FFT SMOOTHING WINDOW SIZE: ! 
READ(*,*) PTS 


OPEN(10, FILE=INFILE, STATUS="OLD" , FORM=" FORMATTED ) 
OPEN (11, FILE=OUTFILE ;STATUS— NEW) 


LASTI=2048 

J=1 

I=0 

PRINT'(/A/)',' Computes ensemble average of turb. inten based on 

local (smoothed) mean:' 

PRINT '(/A/)',' READING INPUT FILE....’ 

PRINT '(A/)',' CAUTION, Watch the number of points in each column 
-.. they should all be about the same length.' 


CONTINUE 
Read the data into array "ENSMBL", each column is a pulse *** 
READ (10,*,END=20,ERR=20) U, PROP 
T=r+1 
But skip the first partial pulse ****** 
IF (J .GT. 1) ENSMBL(I,J-1) =U 
If you have a new pulse, start a new column (prop trigger currently 
set to 1. Volt). *2= 
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IF(I .GE. 20 .AND. PROP .GT. 1.) THEN 
IF (J .EQ. 1) THEN 
WRITE(*,*) 'Skipped first',I,' points (not a full pulse).' 


ELSE 
WRITE(*,*)'No. of points in column',J-1,' is',I 
ENDIF 
kak Find the shortest column *** 


IF(I .LT. LASTI .AND. J .GT. 1) LASTI=I 
IF (J-1 .GE. IPLS) GO TO 28 
J=J3+1 
I=0 
ENDIF 
GO TO 10 
20 CONTINUE 
PRINT '(/A/)',' End of File encountered - deleting last partial 
1 column....' 
J=J-1 
28 CONTINUE 
XN=FLOAT (J-1) 


NPLS=J-1 

c 
WRITE(*,*) 'The final number of pulses is',NPLS,' each of length', 
1LASTI 

Cc 


**x* Perform the FFT smoothing on each column of the array ENSMBL *** 
DO 300 L=1,NPLS 
WRITE(*,%*) 'Smoothing column',L 
DO 310 I=1,LASTI 
UCOL(I) =ENSMBL(I, L) 


310 CONTINUE 
CALL SMOOFT (UCOL, LASTI, PTS) 
kak Write the square error term [{ (U-Usmooth)**2] back into ENSMBL *** 


BO. 320 =), LASTI 
ENSMBL(I,L)=(ENSMBL(I,L) -UCOL(I) ) **2 
320 CONTINUE 
300 CONTINUE 
Cc 
**x* Compute the average local variance and turbulence intensity based on 
* user input value of U - edge (reference velocity), and based on 
* the variance of U from the local smoothed U. Note, the variance is not 
* the true statistical variance. 
DO > 22 I=1,LASTI 
SUM=0. : 
SUMSQ=0. 
DO 30 L=1,J-1 
SUMSQ=SUMSQ+ENSMBL(I,L) 
30 CONTINUE 
VAR (I) =(SUMSQ) / (XN-1.) 
TURBIN (I) =( (SUMSQ/XN) **.5) /UINF 
22 CONTINUE 
zee Write the output ****xx 
DO 40 I=1,LASTI 
XI=FLOAT (TI) 
kak note: variance is currently not output, write array VAR if desired 
WRITE (11,200) XI, TURBIN(T) 
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200 FORMAT (1X,F5.0,F10.5) 
40 CONTINUE 
CLOSE (10) 
CLOSE (11) 
END 
OES SESEZEREREEE REESE SRES ESSER ER ESERRARR ESE SER ERE SE SEER ESE SESE SESE SSE ES ES 


SUBROUTINE SMOOFT(Y,N, PTS) 


Y IS THE DATA ARRAY, N IS THE LENGTH OF THE DATA ARRAY, AND PTS IS 
THE SMOOTHING WINDOW WIDTH 


This subroutine is from the book "Numerical Reciepes" by Press, 
Flannery, Teukolsky and Vetterling, Pub. 1986, Cambridge Univ. Press 
PARAMETER (MMAX=2048) 
DIMENSION Y(MMAX) 


M IS CALCULATED TO BE THE SMALLEST POWER OF 2 WHICH IS GREATER THAN 
THE NUMBER OF POINTS N. THE REMAINDER OF THE ARRAY IS ZERO PADDED. 


NaNaQAN 


M=2 
NMIN=N+2.*PTS 
1 IF(M.LT.NMIN) THEN 
M=2*M 
GO TO 1 
ENDIF 
IF(M.GT.MMAX) PAUSE 'MMAX too small' 
CONST=(PTS/M) **2 
Y1=Y(1) 
YN=Y(N) 
RN1=1./(N-1.) 


REMOVE THE LINEAR TREND 


NaN 


DO 11 J=1,N 
Y (J) =¥ (J) -RN1* (Y1* (N-J)+YN* (J=-1) ) 
a CONTINUE 
IF(N+1.LE.M) THEN 


C 
C ZERO PAD 
e 
DO 12 J=N+1,M 
Y(J)=0. 

12 CONTINUE 

ENDIF 
e 
C TRANSFORM TO FREQUENCY DOMAIN 
Cc 

MO2=M/2 

CALL REALFT(Y,MO2,1) 

Y(1)=¥ (1) /MO2 

FAC=1. 
@ 
é WRITE(11,*) ' BIN (I) Y(F)! 
Ce DO 99 I=1,M 
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WRITE( 11.2%) 1,¥ (1) 
99 CONTINUE 


WINDOW FUNCTION - (USES THE PARABOLIC WELCH WINDOW) 


ANANAAAN 


WRITE(11,*) ' WINDOW FUNCTION: ! 
DO’ 13 J=1,M02—1 
K=2*J+1 
IF (FAC.NE.0.) THEN 
FAC=AMAX1(0., (1.-CONST*J**2) /MO2) 
Cc WRITE (145%) K,K+1, XC 
Y (K) =FAC*Y (K) 
Y¥ (K+1) =FAC*®Y (K+1) 
ELSE 
Y (K)=0. 
Y (K+1) =0. 
ENDIF 
13 CONTINUE 
Cc 
C HANDLES THE LAST POINT 
Cc 
FAC=AMAX1(0., (1.7-0.25*PTS**2) /MO2) 
Y(2)=FAC*Y (2) 
WRITE(11,%*) ‘ WINDOW x DATA:' 
DO 98 I=1,M 
Weibel, =) 1,Y¥(L) 
98 CONTINUE 


TRANSFORM BACK TO TIME DOMAIN 
CALL REALFT(Y,MO2,-1) 


RESTORE LINEAR TREND 


NymaAaaQq aNANAANANANAN 


DO 14 J=1,N 
¥ (J) =RN1* (Y1* (N-J) +YN* (J=-1) )+Y(J) 

CONTINUE 

RETURN 

END 


CREEK KRKERKEKKEKKKKEKEKREKERKEREKKKEKEKKKEERKKRRREKEKKERREKRERRRRKRERERKEKREKE 


SUBROUTINE REALFT (DATA,N, ISIGN) 
CALCULATES THE FOURIER TRANSFORM OF A SET OF 2*N REAL VALUED DATA 
POINTS (ARRAY 'DATA'). ISIGN IS 1 FOR THE FOURIER TRANSFORM (TIME 
TO FREQ. DOMAIN) AND ISIGN IS -1 FOR INVERSE TRANSFORMATIOM. 
THIS SUBROUTINE SETS UP THE DATA TO BE EFFICIENTLY TRANSFORMED BY 
SUBROUTINE 'FOUR1'. 


~~ 
> 


This subroutine is from the book "Numerical Reciepes" by Press, 
Flannery, Teukolsky and Vetterling, Pub. 1986, Cambridge Univ. Press 
REAL*8 WR,WI,WPR,WPI,WTEMP, THETA 
DIMENSION DATA(*) 
THETA=6. 28318530717959D0/2.0D0/DBLE (N) 
C1=0.5 


NAAANAANAANANANANA 
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FORWARD TRANSFORM BY FOUR1 


NaaQN0 


IF (ISIGN.EQ.1) THEN 
C2=-0.5 
CALL FOUR1(DATA,N, +1) 
ELSE 


SET UP FOR INVERSE TRANSFORM 


aqgaaqgn 


C2=0.5 
THETA=-THETA 

ENDIF 

WPR=-2.0D0*DSIN(0.5DO*THETA) **2 

WPI=DSIN (THETA) 

WR=1.0D0+WPR 

WI=WPI 

N2P3=2*N+3 

DO 11 I=2,N/2+1 
I1=2*I-1 
I2=I1+1 
I3=N2P3-I2 
14=I3+1 
WRS=SNGL (WR) 
WIS=SNGL (WI) 
H1R=C1* (DATA(I1)+DATA(TI3) ) 
H1I=C1* (DATA (I2) -DATA(I4) ) 
H2R=-C2* (DATA (I2)+DATA(I4)) 
H2I=C2* (DATA(I1) -DATA(I3) ) 
DATA (1I1) =H1R+WRS *H2R-WIS*H2I 
DATA (I2) =H1I+WRS*H2I+WIS*H2R 
DATA (13) =H1R-WRS *H2R+WIS*H2I 
DATA (14) =-H1I+WRS*H2I+WIS*H2R 
WIEMP=WR 
WR=WR*WPR-WI *WPI+WR 
WI=WI*WPR+WTEMP*WPI+WI 

ala CONTINUE 

IF (ISIGN.EQ.1) THEN 
H1R=DATA (1) 
DATA (1) =H1R+DATA (2) 
DATA (2) =H1R-DATA (2) 

ELSE 
H1R=DATA(1) 
DATA (1)=C1*(H1R+DATA(2)) 
DATA (2)=C1* (H1R-DATA(2)) 


Cc 
C INVERSE TRANSFORM 
Cc 
CALL FOUR1(DATA,N,-1) 
ENDIF 
RETURN 


END 


CRRA RKRAKKKKKKKKKKKKKEKRKEKRKEKKKKKEKKKKKKKKKKkkKkkKKKKKKKKKKKKKKXKKEKE 


SUBROUTINE FOUR1(DATA,NN, ISIGN) 
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COMPUTES THE COMPLEX DFT OR INVERSE DFT OF ARRAY ‘DATA‘' OF 
BENGIH NN’. ISIGN IS 1 FOR) THE DFY AND -1 FOR THE INVERSE 
DET. 


This subroutine is from the book "Numerical Reciepes" by Press, 


REAL*8 WR,WI,WPR,WPI,WTEMP, THETA 
DIMENSION DATA(*) 
N=2*NN 
J=1 
DO 11 I=1,N,2 
IF(J.GT.1I) THEN 
TEMPR=DATA (J) 
TEMPI=DATA (J+1) 
DATA (J) =DATA(TI) 
DATA (J+1) =DATA(I+1) 
DATA (I) =TEMPR 
DATA (I+1) =TEMPI 
ENDIF 
M=N/2 
1 IF ((M.GE.2).AND.(J.GT.M)) THEN 
J=JI-M 
M=M/2 
GO TO 1 
ENDIF 
J=I+M 
11 CONTINUE 
MMAX=2 
2 IF (N.GT.MMAX) THEN 
ISTEP=2 *MMAX 
THETA=6. 28318530717959D0/ (ISIGN*MMAX) 
WPR=-2.D0*DSIN(0.5DO*THETA) **2 
WPI=DSIN (THETA) 
WR=1.D0 
WI=0.DO 
DO 13 M=1,MMAX,2 
DO 12 I=M,N,ISTEP 
J=I+MMAX 
TEMPR=SNGL (WR) *DATA(J) -SNGL(WI) *DATA (J+1) 
TEMPI=SNGL(WR) *DATA(J+1)+SNGL(WI) *DATA(J) 
DATA(J)=DATA(I)-TEMPR 
DATA (J+1)=DATA (I+1)-TEMPI 
DATA (I) =DATA(I)+TEMPR 
DATA (I+1)=DATA (I+1)+TEMPI 
2 CONTINUE 
WTEMP=WR 
WR=WR*WPR-W1I*WPI+WR 
WI=W1*WPR+WTEMP*WPI+WI 
13 CONTINUE 
MMAX=ISTEP 
GO TO 2 
ENDIF 
RETURN 
END 
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Flannery, Teukolsky and Vetterling, Pub. 1986, Cambridge Univ. Press 
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USER’S GUIDE FOR PROGRAM PSD 


1. Input Instructions 


Insure that the executable file for program PSD (PSD.EXE) 
and the input data file are in the same subdirectory of 
a disk drive with sufficient space remaining to write the 
output file. Note that even though program CORREX does 
not work with an ensemble, the program is currently set 
up to read a data file with two inputs such as U and PROP 
described in section B.1. The PROP variable is not used 
and is only a "place holder" in the read statement. The 
read statement is easily modified to suit the user's 
needs. 


Initiate program PSD (type "PSD" at the DOS prompt). 


The program will prompt the user for the name of the 
input data file. Enter the full name of the input data 
file including the extension. 


The program will prompt the user for the name of the 
output file which will be created by the program. Enter 
the desired file name. Do not use the name of an 
existing file or the program will terminate. 


The program will prompt the user for the sample rate of 
the date. Enter the sample rate of the data contained 
in the input file. 


The program will prompt the user to select one of the 
seven tapering window options from Table IV-4. Enter the 
number of the desired option. (Windows 5 and 6 are good 
compromises. ) 


The program will prompt the user to choose whether or not 
to remove the stationary mean prior to computing the 
spectral estimate. Enter "Y" for yes (remove the mean) 
or "N" for no (do not remove the mean). If the data has 
a large mean value (relative to the magnitude of the 
oscillations) then it is recommended that the mean be 
removed. 


The program will prompt the user to enter the desired 
number of frequency bins -(M) and half of the number of 
overlapping segments (K). Enter the selected values of 
M and K separated by a comma. WARNING: The values for 
M and K may NOT be selected arbitrarily! The allowable 
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values of M and K are dependent on each other and on the 
total number of data points in the time series (N) 
according te the equation N = (2K + 1)M where M must be 
an integer power of two and N is equal to or slightly 
less than an integer power of two. The tradeoff is too 
obtain an adequate frequency resolution (large M) and 
still minimize the variance of the spectral estimate 
sufficiently (large K). Table IV-3 presents some typical 
values of M and K for different record lengths N. 
Usually, the total number of data points (N) will not be 
an integer power of two. The user may choose to read a 
lesser number of data points (reasonable if N is large, 
say greater than 4096). Or the user can add zeros to the 
end of the data file to achieve the next larger integer 
power of two value for N. Note that the actual number 
of points read by the program is determined by the user 
input values of M and K. 


- The program will then compute the spectral estimate and 
write the output file. 


2. Output Description 

The output file consists of five columns. The first 
column is the frequency. The second column contains the 
magnitude of the power at each listed frequency. The third 
column presents the relative power in decibels (dB) computed 
such that the peak power is zero QB. The fourth column 
presents a running total of the power (magnitude) and the 
fifth column lists the percent running total power. The 


following pages present the source code for program ENSEMBL. 
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PROGRAM PSD 


REMARK KEKE ERERAKAERKEKEREKRAEKKAKEKRKEKKKRKEKKKEKAKKKAEKKEKKEKKKKARKEESE 
This program computes an estimate of the PSD of a real discrete time 
series using the classical averaged periodogram method with 50% overlapped 
segments for (nearly) minimum variance. The main program is a driver to 
set up the input for the call to SUBROUTINE SPCTRM which does most of 

the work including reading the data. 


Input file is a sequential file containing values of velocity and PROP. 
(Such as the output of the Labtech Notebook A/D conversion software. Time 
is implied by the sample rate. 


The key variables are defined below: 


M = Number of freq. bins and half the number of points per segment. 
Large M implies high resolution in frequency. (MUST BE A POWER OF 2) 
K = 1/2 of the number of segments. Large K implies low PSD variance. 
K is given by: Kaze NT ( (HM 172) 
N = Total number of data points (power of 2). N = (2K+1)M 


TYPICAL EXAMPLES: N M K 
al2 32 z, 
64 3 

128 1 

1024 32 25 
64 7 

P28 3 

2048 64 iS 
12s 7 

25:6 3 

4096 64 31 
128 pls 

256 U 

S292 64 63 
128 a 

256 ale) 

Dil u 


SMPLRT = Sample rate of the input data. 

FREQ = Frequency value in HZ of a particular frequenct bin. 

DFREQ = The width of each frequency bin. 

Wl and W2 = Workspace arrays for reading the input data sequence from 
unit 9: 

P = Output PSD array in (ft/sec)**2 per Hz 

DBP = Output PSD array in decibels per Hz 

Tl and T2 = Dummy read variables. 

WINDOW = The current FFT window function value (see function subroutine). 


NQNANANAANANQNANANANANANANANANANAANANNANANANANANNANAANANANAANANANANNANANANANARAANANRANAANANARA 


WOPT = Window selection parameter: WOPT WINDOW 
1 Parzen (triangle) 
2 rectangular (none) 
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Welch 

Van Hann 
Hamming 

3 term B-H 
4 term B-H 


SAU & W 


Cc 

c 

Cc 

cS 

Cc 

c 

C Note on dimensions: P and PSD are currently dimensioned to ((2**L) + 1) 
C frequency bins where L=9. (That's 513, so the MAXIMUM frequency 

C resolution is currently 513 bins.) W2 and Wl are dimensioned to 

C (2**L) and 4*(2**L) respectively. (That's 512 and 2048.) If you need 

C more resolution, the next higher dimensions are 1025 for P and DBP, and 
C 1024 and 4096 for W2 and Wl respectively. 

c 
CG 
€ 
6 


Programmed by: Donald K. Johnson 
August 1988, NPS Monterey, CA 
BREAKER AAAAAAEAERAKAEAEKEAERRAEEAAKKAAEREAAKKAKEAAEKREAREKEKEKKKERAEARRARARKKKK 
REAL SMPLRT, DFREQ, FREQ, T2, Tl 
INTEGER 2,0 ,K,M,WOPT 
DIMENSION P(513) ,W1(2048) ,W2 (512) ,DBP(513) , TPWR(513) , PCTPWR(513) 
LOGICAL OVRLAP 
CHARACTER*12 INFILE,OUTFILE 
CHARACTER*1 RM 


Cc 
PRINT '(/A/)',' Welch Periodogram Spectral Estimation Program' 
c 
PRINT '(/A\)',' ENTER THE INPUT DATA FILE NAME: ! 
READ (*,'(A)') INFILE 
Cc 
PRINT "(7A\)",' ENTER THE OUTPUT FILE NAME: ! 
READ (*,°(A)*) OUTFILE 
Cc 
PRINT '(/A\)',' ENTER THE SAMPLE RATE OF THE DATA: ! 
READ (*,*) SMPLRT 
Cc 
WRITE(*,'(/A\)')' ENTER FFT "WINDOW" OPTION (1 thru 7): ' 
READ(*,*) WOPT 
Cc 
PRINT '(/A\)',' DO YOU WANT TO REMOVE THE MEAN FIRST (Y OR N)? ! 
READ (*, (A) ) RM 
(Gs 
PRINT “C/AN) > ENTER THE NO. OF FREQ BINS (M) AND HALF OF THE NO. 
1 OF SEGMENTS (K): ! 
READ (*,*)'M,K 
M4=4*M 
Cc 
OPEN (UNIT=9, FILE=INFILE, STATUS='OLD') 
OPEN (UNIT=11, FILE=OUTFILE, STATUS='NEwW') 
Cc 
zx*x if selected, remove the mean and write the u' into scratch file unit 10 


KUNIT=9 
IF (RM .EQ. 'Y' .OR. RM .EQ. 'y') THEN 
OPEN (UNIT=10, STATUS='SCRATCH') 

NPTS=(2*K+1) *M 
PRINT '(/A\,I5)',' NPTS = ',NPTS 
SUMD1=0. 
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compute the sum 
DO 2 I=1,NPTS 

READ(9,*) U 

READ(9,*) U, PROP 

SUMD1=SUMD1+U 
CONTINUE 
REWIND(9) 
compute the mean 
AVG1=SUMD1/NPTS 
PRINT '(/A\,F12.4)',' The removed mean is: ‘,AVG1 
compute u' (remove the mean) 
DO 4 I=1,NPTS 

READ(9,*) U 
READ(9,*) U, PROP 

UPRM=U-AVG1 

WRITE(10,*) UPRM, PROP 
CONTINUE 
REWIND(9) 
REWIND(10) 
change the read file unit number from 9 to 10 (passed to subroutine) 
KUNIT=10 

ENDIF 


DFREQ=SMPLRT/ (2*M) 
OVRLAP=.TRUE. 


The subroutine itself contains the actual read statements 


PRINT ‘'(/A/)',' Estimating spectrums. .‘ 
CALL SPCTRM(P,M,K,OVRLAP,WOPT,KUNIT,W1,W2,T1,12,70BP, TPWR, Polen) 
WRITE(*,*) ‘SPECTRUM OF °, INFILE,’ IS IN PELE OU re ree 
DO 11 J=1,Mt+1 
FREQ=J*DFREQ ~- DFREQ 
WRITE(11, '(1X,F10.1,F17.3,F10.3,Fi7.3,F8.2)") FREG eee eee 


1 TPWR(J),PCTPWR(J) 


CONTINUE 
CLOSE(9) 

CLOSE(10) 

CLOSE(11) 

PRINT '(/,/A)',' ALL DONE! ! 
END 


KR REA K KK HRKREKEKRK KR KRRARRRAEKRRKERKREKRKAEREKAKREKEKKAREKRKKKKAEKAEKRKRKRKEKKRKEKKEKKEKEKR 


AQAANANANAAANAANAANA 


SUBROUTINE SPCTRM(P,M,K,OVRLAP,WOPT, KUNIT,W1,W2,T1,T2,DBP,TPWR, 


+PCTPWR) 


This subroutine represents a modified version of a spectral analysis 
program published in "Numerical Recipes" by Press, Flannery, Teukolsky, 
and Vetterling, pub 1986, Cambridge Univ. Press. 


Reads data from unit 9 and returns as P(J) the data's power (mean square 
amplitude) at frequency (J-1)/(2*M) which is then converted to HZ in the 
Iain program, for J=1 to M based on (2*K+1)*M data points. 


modified by: Donald K. Johnson 
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August 1988, NPS Monterey, CA 


LOGICAL OVRLAP 
INTEGER WOPT 
DIMENSION P(M+1),W1(4*M) ,W2(M) ,DBP(M+1),TPWR(M+1) , PCTPWR(M+1) 
MM=M+M 
M4=MM+MM 
M44=M4+4 
M4 3=M4+3 
DEN=0. 
FACM=M-0.5 
FACP=1./(M+0.5) 
PI2=8.*ATAN(1.) 
SUMW=0. 
compute the sum of the squared window function for later normalizing *** 
the window functions are given in the function subroutine below. 
DO ii) J=1,MéMM 

SUMW=SUMW+WINDOW (J ,MM,WOPT) **2 
CONTINUE 
initialize the spectrum to zero *** 
DO 12 J=1,M 

P(J)=0. 
CONTINUE 
initialize the "save" half buffer *** 
IF (OVRLAP) THEN 


for overlapping, W2 is the velocity and, T2 is a dummy 
READ (KUNIT,*) (W2(J),T2,J=1,M) 
READ (KUNIT,*) (W2(J),J=1,M) 


ENDIF 
loop over the data segments in groups of two. Get two complete 
segments into the workspace. *** 
DO 18 KK=1,K 
DOT 1 SeJORE——1, 0, 1 
IF (OVRLAP) THEN 
DO 13 J=1,M 
W1 (JOFF+J+J) =W2 (J) 
CONTINUE 


for overlapping, W2 is the velocity and T2 is a dummy 
READ (KUNIT,*) (W2(J),T2,J=1,M) 
READ (KUNIT,*) (W2(J),J=1,M) 


JOFFN=JOFF+MM 
DO 14 J=1,M 

W1 (JOFFN+J+J) =W2 (J) 
CONTINUE 

ELSE 
for no overlap, Wl is the velocity and Tl is a dumny 

READ (KUNIT,*) (W1(J),T1,J=JOFF+2,M4,2) 

READ (KUNIT,*) (W1(J),J=JOFF+2,M4,2) 


ENDIF 
CONTINUE 
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C *** apply the window function defined above to the data *** 
DO 16 J=1,MM 
J2=J+I 
W=WINDOW(J,MM, WOPT) 
W1(J2)=W1(J2) *W 
W1(J2-1)=W1(J2-1) *W 


16 CONTINUE 
C *** Fourier transform the windowed data, then sum results into previous 
@ segment *** 


CALL FOUR1(W1,MM,1) 
P(1)=P(1)+W1(1) **2+W1(2) **2 
DO 17 J=2,M 
J2=I+J 
P(J)=P(J)+W1(J2) **2+W1(J2-1) **2 
+W1(M44-J2) ##2+W1(M43-J2) **2 
17 CONTINUE 
P(M+1)=P(M+1)+W1(MM+1) **2+W1 (MM+2) #*2 
DEN=DEN+SUMW 


18 CONTINUE 
C *** correct normalization *** 
DEN=M4 *DEN 


C *** normalize the output *** 
DO 19 J=1,M+1 
P(J)=P(J)/DEN 


1G.) CONTINUE 
C*** compute the relative PSD in decibels and the running total power *** 
PMAX=P(1) 


DO 50 J=1,M 
IF(P(Jt+1) .GT. PMAX) PMAX=P(J+1) 
S50 CONTINUE 
SUMPWR=0. 
DO 60 J=1,M+1 
IFO (PP (0) .EQ. 0.) THEN 
WRITE(*,*) "WARNING: P(J) has a zero value. P(J) set to -999 
1999.0" —4t ad cf!,J,° -.- CONEDTMUING 
DBP (J) =-999999. 
ELSE 
DBP(J)=20. *ALOG10(P(J) /PMAX) 
SUMPWR=SUMPWR+P(J) 
TPWR(J)=SUMPWR 
ENDIE 
60 CONTINUE 
**x* compute the percentage of total power at freq. J relative to the 
* total power summed over all freqs. 
DO 70 J=1,M+1 
PCTPWR(J)=100. *TPWR(J) /TPWR(M+1) 
70 CONTINUE 
RETURN 
END 


CARRERE EAEKEEEEKEKEEEKER 


FUNCTION WINDOW(J,MM,WOPT) 


These window functions come from the paper "On the Use of Windows 
for Harmonic Analysis with the Discrete Foutier Transform", by 
Fredric J. Harris, IEEE Proc., Vol. 66, NO-.1-9Jan i970- 


+e + © 
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J = The current time index. 
MM = The total number of points to be windowed. 
WOPT = The selected window option. 


» + 4+ & 


INTEGER WOPT 
PI2=8.*ATAN(1.) 
M=MM/2 
FACM=M-0.5 
FACP=1./(M+0.5) 
IF(WOPT .EQ. 1) THEN 
Cc -- PARZEN (TRIANGULAR) WINDOW: 
WINDOW=(1.-ABS (((J-1) -FACM) *FACP) ) 


* 


(ay 
ELSEIF(WOPT .EQ. 2) THEN 

Cc -- RECTANGULAR WINDOW: ie. no window 
WINDOW=1. 

Cc 


ELSEIF(WOPT .EQ. 3) THEN 
C -- WELCH (PARABOLIC) WINDOW: 
WINDOW=(1.-(((J-1)-FACM) *FACP) **2) 


ELSEIF(WOPT .EQ. 4) THEN 
Cc -- VAN HANN (SQUARED COSINE) WINDOW: 
WINDOW=0 .5-0.5*COS ( PI2* (FLOAT (J-1) /FLOAT (MM) )) 


ELSEIF(WOPT .EQ. 5) THEN 
e -- HAMMING (RAISED COSINE) WINDOW: 
WINDOW=0.538-0.462*COS (PI2* (FLOAT (J-1) / FLOAT (MM)) ) 


Cc 
ELSEIF(WOPT .EQ. 6) THEN 
6 -- BLACKMAN-HARRIS 3 TERM MINIMUM WINDOW: 
WINDOW =.42323-.49755*COS (PI2* (FLOAT (J-1) /FLOAT (MM) ))+ 
1 .07922*COS (PI2*2. * (FLOAT (J-1) /FLOAT (MM) )) 
Cc 
ELSEIF(WOPT .EQ. 7) THEN 
Cc -- BLACKMAN-HARRIS 4 TERM MINIMUM WINDOW: 
WINDOW =.35875-.48829*COS (PI2* (FLOAT (J-1) /FLOAT(MM)))+ 
i. .14128*COS (PI2*2.* (FLOAT (J-1) /FLOAT(MM)))- 
2 .01168*COS (PI2*3.* (FLOAT (J-1) /FLOAT (MM) )) 
Cc 
ELSE 
WRITE(*,*) 'INVALID WINDOW SELECTION, TERMINATE...' 
STOP 
ENDIF 
RETURN 
END 


Ce RRR ERE RK RKERRRRERR RE REKKERREEAKARRERERREAR 
SUBROUTINE FOUR1(DATA,NN,ISIGN) 


Subroutine from "Numerical Recipes" by Press, Flannery, Teukolsky, 
and Vetterling, pub 1986, Cambridge Univ. Press. 


AYANAAN 


REAL*8 WR,WI,WPR,WPI,WTEMP, THETA 
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DIMENSION DATA(*) 
N=2*NN 
J=1 
DO 2b r= NN 
IF(J.GT.1) THEN 
TEMPR=DATA (J) 
TEMPI=DATA (J+1) 
DATA (J) =DATA (TI) 
DATA (J+1)=DATA(I+1) 
DATA(I)=TEMPR 
DATA (I+1)=TEMPI 
ENDIF 
M=N/2 
TF ((M.GE-2).AND<(J .GI.M))) tHe 
J=J-M 
M=M/2 
GO TO l 
ENDIF 
J=I+M 
CONTINUE 
MMAX=2 
IF (N.GT.MMAX) THEN 
ISTEP=2 *MMAX 
THETA=6.28318530717959D0/ (ISIGN*MMAX) 
WPR=-2.D0*DSIN(0.5DO*THETA) **2 
WPI=DSIN (THETA) 
WR=1.D0 
WI=0.D0 
DO 13 M=1,MMAX,2 
DO 12 I=M,N,ISTEP 
J=I+MMAX 
TEMPR=SNGL(WR) *DATA(J) -SNGL(WI) *DATA(J+1) 
TEMPI=SNGL(WR) *DATA(J+1)+SNGL(WI) * DATA (J) 
DATA (J) =DATA(I) -TEMPR 
DATA (J+1)=DATA(I+1) -TEMPI 
DATA (I) =DATA(I)+TEMPR 
DATA (I+1)=DATA(I+1)+TEMPI 
CONTINUE 
WTEMP=WR 
WR=WR*WPR-WI*WPI+WR 
WI=WI*WPR+WTEMP*WPI+WI 
CONTINUE 
MMAX=ISTEP 
GO 70 22 
ENDIF 
RETURN 
END 
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USER’S GUIDE FOR PROGRAM ENSPSDAV 


1. Input Instructions 


Insure that the executable file for program ENSPSDAV 
(ENSPSDAV.EXE) and the input data file are in the same 
subdirectory of a disk drive with sufficient space 
remaining to write the output file. 


Initiate program ENSPSDAV (type "ENSPSDAV' at the DOS 
prompt). 


The program will prompt the user for the name of the 
input data file. Enter the full name of the input data 
file including the extension. 


The program will prompt the user for the name of the 
output file which will be created by the program. Enter 
the desired file name. Do not use the name of an 
existing file or the program will terminate. 


existing file or the program will terminate. 


The program will prompt the user for the sample rate of 
the data. Enter the sample rate of the data contained 
in the input file. 


The program will then read the input data file. Remember 
to observe the column lengths of the array ENSMBL to 
insure that the ensemble is being constructed properly. 


Next, the program will prompt the user to choose whether 
or not to remove the mean of each individual segment 
prior to computing the spectral estimate of the segment. 
Enter "Y" for yes (remove the mean) or "N" for no (do not 
remove the mean). If the data has a large mean value 
(relative to the magnitude of the oscillations) then it 
is recommended that the mean be removed. 


The program will prompt the user to select one of the 
seven tapering window options from Table IV-4. Enter the 
number of the desired option. (Windows 5 and 6 are good 
compromises). 


Then the program will prompt the user to select one of 
three options for segmenting the data. Enter the number 
of the selected option. Option 12 allows the user to 
enter the number of samples between the beginning of each 
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segment. Option 23 allows the user to select the number 
of segments which span the ensemble (i.e., the number of 
spectral estimates across the ensemble of pulses). 
Option 3 allows the user to enter the sample number of 


the mid-point location of a single segment (i.e., the 
location in time of a single spectral estimate within the 
ensemble). All segments are 256 samples wide. 


- The program will then compute the spectral estimates of 
the selected segments and write the output file. 


¢ The program will prompt the user to choose whether to end 
operation or continue working with this same ensemble of 
data. If the user chooses to continue, the program 
prompts for another output file name and the process 
begins again. 

2. Output Description 
The output file consists of some heading information 
and five columns of spectral data. The heading information 


summarizes the options selected (segmenting option, tapering 


window option, the mean removal option and the input data file 


name). Then the location of the first segment spectral 
estimate is identified. The first column of data is the 
frequency. The second column contains the magnitude of the 


power at each listed frequency. The third column presents the 
relative power in decibels (dB) computed such that the peak 
power is zero @dB. The fourth column presents a running total 
of the power (magnitude) and the fifth column lists the 
percent running total power. Then the location of the next 
segment is identified and the five columns of spectral data 
are repeated for this next segment, repeating until the 


spectral results for each ensemble segment have been reported. 
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The following pages present the source code for program 


ENSEMBL. 
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$DEBUG 


PROGRAM PSDAVG 


RRMA KRRERAKAKAAEEREREKRAEEERERKRAAERKRKEKAKEEKERAKAKKEKKKEKAKKKKKKKKKKKKK 


+e eee ee eee ee ee OU Uh ll UU ll lt ll Oe OO OO ee Oe Oe OOO 


Program to perform spectral analysis on short time segments of an 
ensemble of non-stationary turbulent boundary layer velocity data to 
obtain the time varying spectral characteristics of the data. 


Input data file is assumed to be in the format of the A/D conversion 
system, ie. a sequential file with alternating values of velocity (U) 
and propellor pulse (PROP), starting with velocity. 


The user is given three options on how to segment the ensemble, 

by specifying the no. of points between the beginning of each segment, 
the no. of segments across the ensemble, or the center point of a 
single segment. 


After segmenting per the users choice, and removing the mean (if 
selected), the program builds array Wl which contains M consecutive 

time points from a single pulse. This time slice is windowed using the 
spectral window selected.The program, then computes a PSD estimate on the 
time slice. This process is repeated for the corresponding time slice 

in each propeller pulse. Then the PSD of each pulse time slice is 
averaged to get the final PSD for the segment. This PSD is output 

and then the entire sequence is repeated on the next segment. 

The next segment is shifted IWIDE time points, but is always M points 
long. 


Key variables are defined below: 


ENSMBL = The primary data array which contains up to 50 pulses (columns) 
each of length up to 2048 points. The pulses are synchronized in 
time (an ensemble) such that each point in a pulse "lines up" 
with the others. 

U = The measured velocity as read from the input file into array ENSMBL. 

PROP = The propeller pulse timing marker signal. Also read from the input. 

IPLS = Number of pulses to be read from the input file. 

LASTI = Number of points in the columns of ENSMBL. (lie. rows in ENSMBL) 


IWIDE Number of time points from the beginning of the current segment 
to the beginning of the next segment. 
NSEG Number of time slice segments in the ensemble. 


NPLS = The actual total number of pulses (columns in ENSMBL) after the 
input data has been split into an ensemble. 
ISTART = Starting point number of the segment. 
PSD Related Variables: 


M = Number of points per segment (MUST BE A POWER OF 2). Also the number 
of points in each FFT. (Currently fixed at 256 points) 


WINDOW = The current FFT window function value (see function subroutine). 


WOPT = Window selection parameter: WOPT WINDOW 
al Parzen 
2 Rectangular 
3 Welch (parabolic) 
4 von Hann 
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4 term B-H 


SAWN 


See the following reference for details on the window functions: 
"On the Use of Windows for Harmonic Analysis with the DFT" by 
Fredric J. Harris, IEEE Proc., Vol. 66, No.1, January 1978 


SMPLRT = Sample rate of the input data. 

FREQ = Frequency value in HZ of a particular frequenct bin. 

DFREQ = The width of each frequency bin. 

Wl = A complex workspace array for the FFT subroutine FOUR]. 

P = Output PSD array in (ft/sec)**2 per Hz 

DBP = Output PSD array in decibels per Hz 

TPWR = The running total power summed from zero frequency to the Jth freq. 

PCTPWR = Percentage of total power summed from zero frequency to the Jth 
freq. relative to the total power summed across all freqs. 


Note on dimensions: P is currently dimensioned to M/2+1 (129) frequency 
bins. Wi is dimensioned to 2*M because it is a complex array. The even 
elements contain the imaginary part (all zero in this case) and the 

odd elements contain the real part (the velocity data). 


Programmed by: Donald K. Johnson 
August 1988, NPS 


REE KEKEAEEKEKREKKEKKEKEKREKKEKKKKEKeEeKKKKKEKKEKEKEKEEKKEeEKKEKEEKEE KEKE KEKKEKKEKEESE 


PARAMETER (MMM=256 , MMMO2=128) 
DIMENSION ENSMBL(2048,50) ,P(MMMO2+1) ,W1(2*MMM) , DBP (MMMO2+1) 
DIMENSION TPWR(MMMO2+1) , PCTPWR(MMMO2+1) 

INTEGER I, IPLS,NPLS,J,L,LL,NPTS, IWIDE, WOPT 

REAL FREQ, DFREQ, SMPLRT 

CHARACTER*12 INFILE, OUTFILE 

CHARACTER*1 RM,CONT 


PRINT '(/A/)',' Spectral Analysis of Pulse Segments' 
PRINT '(/A\)',' ENTER THE INPUT DATA FILE NAME: ! 

READ (*,'(A)') INFILE 

PRINT '(/A\)',*' ENTER THE OUTPUT DATA FILE NAME: ' 

READ (*,°*(A)") OUTFILE 

PRINT '(/A\)',' ENTER THE NUMBER OF PROP PULSES YOU WANT: '! 
READ (*,*) IPLS 

IPLS=50 

PRINT ‘(/A\)','* ENTER THE SAMPLE RATE OF THE DATA: '! 


READ(*,*) SMPLRT 


OPEN (10, FILE=INFILE, STATUS='OLD') 
OPEN (11, FILE=OUTFILE, STATUS='NEW') 


LASTI=2048 
J=1 
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I=0 
PRINT'(/A/)',' Reading input file...Watch the no. of points dneeae 
lh column! * 


CONTINUE 
Read the data into array “ENSMBL", each column is a pulse *** 
If an EOF is encountered then exit the loop and jump to 20. If no 
EOF is encountered then jump to 28 
READ (10,*,END=20,ERR=20) U,PROP 
I=I+1 
But skip the first partial pulse ***##* 
IF (J .GT. 1) ENSMBL(I,J-1)=U 


If you have a new pulse, start a new column (pulse trigger currently 


set to 1 volt. 
IF(I .GE. 20 .AND. PROP <GT. 22)” THEN 
PP (9 =. BO. 1) THEN 
WRITE (*,*) 'Skipped first',I,' points (not a full pulseyee 
ELSE 
WRITE(*,%*) 'No. Of points in column’ )J=1)) sist 
ENDIF 
Find the shortest column *** 
IF(I .LT;: LASTI .«~AND. J <GT. 1S LAST I“: 
Quit loop if you've read enough pulses *** 
IF (J-1..GE.2PES) GO Te 25 
J=J+1 
I=0 
ENDIF 
Curse me, a go to statement! I never took structured programming. *** 
GO TO 10 
CONTINUE 
Since an EOF was encountered, write error message and decrement J to 
ignore the last partial pulse. This is the normal path. 
WRITE(*,'(/A/)')' End of File encountered - deleting last parti 
lal column....' 


J=J-1 
CONTINUE 
XN=FLOAT (J-1) 
NPLS=J-1 
WRITE(*,'(A,I3,A,16)')' The final number of pulses is',NPLS,' , ea 


leh of lengtn *;,GASTI 


WRITE(*,'(/,/A\)')' DO YOU WANT TO REMOVE THE MEAN BEFORE COMPUTIN 
1G PSDs? {Yor N)z 
READ(*,'(A)') RM 


WRITE(*,'(/A\)')' ENTER FFT "WINDOW" OPTION (1 thru 7): ' 
READ(*,*) WOPT 


Menu of options for segmenting ensembles *** 
WRITE(*,'(/A/)')' How do you want to segment the ensemble of pulse 


isi? 


WRITE(*,'(A)')' OPTION DESCRIPTION'! 
WRITE(*,'(A)')' i Select the time width between each PSD' 
WRITE(*, °(A) *)' 2 Select the total number of PSDs (segmen 
its)" 
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WRITE(*,'(A)')! 3 Select the center time point for a sing 


ive Psp! 
WRITE (*,'(/A\)')' ENTER THE OPTION NUMBER OF CHOICE: ' 


READ(*,*) IOPT 


z*x* Note that all segment lengths (M) are 256 points wide. Smaller width 


* separations (IWIDE) cause adjacent segments to be overlapped. 
x 
zz option 1 wx 

TF (OPT .EQ. 1) THEN 


WRITE(*,'(/A\)')' ENTER THE TIME WIDTH (IN NO. OF POINTS) BETW 


1LEEN EACH PSD: ' 
READ(*,*) IWIDE 


aah this scheme will probably cause you to pick up some zeros after 

* the end of the data - so the last segment may have a discontinuity 
NSEG=(LASTI-M)/IWIDE + 1 
ISTART=1 
M=256 


zeke option 2 *** 
EFLSEIF(LOPYT .EQ. 2) THEN 


WRITE(*,'(/A\)')' ENTER THE TOTAL NUMBER OF PSDs (SEGMENTS) AC 


1ROSS THE ENSEMBLE: ' 
READ(*,*) NSEG 


hae you will probably miss a few points at the end of the pulse 
IWIDE=(LASTI-M) / (NSEG-1) 
ISTART=1 
M=256 


* 


zeke option 3 *** 
ELSETE (IOPT <EQ. 3) THEN 


WRITE(*,'(/A\)')' ENTER THE TIME POINT NUMBER (CENTER POINT) 0 
1F A SINGLE 256 POINT TIME SLICE (at least 128 points from eit 


2her end of the pulse): ' 
READ(*,*) ITIME 
ISTART=ITIME-127 
IWIDE=256 
NSEG=1 
M=256 
EESE 
WRITE(*,'(/A/)')' OOPS! - INVALID OPTION CHOICE, TRY AGAIN' 
GO TO 25 
ENDIF ; 
IF(IWIDE .GT. M) THEN 
WRITE(*,'(/A/)')' SEGMENT WIDTH .GT. 256 (TOO LARGE).' 
GO, TO 25 


ENDIF 
* 


*x*x Compute the spectral estimate of an ensemble of pulses in data segments 


* that are M points wide and IWIDE points apart. 
* First initialize some constants. 

ZERO=0.0 

MO2=M/2 

DFREQ=SMPLRT/M 

MM=M+M 
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M44=MM+4 
M43=MM+3 
SUMW=0. 
sum of square of window function, used for normalizing psd later *** 
DO 45 J=1,M 
SUMW=SUMW+WINDOW(J,M, WOPT) **2 
CONTINUE 


write some heading information to the begining of the output file *** 
WRITE(11,'(A,I3)')' Selected options: Segment option: ‘,IOPT 


WRITE(11,*°(A,13) "92 Window option: ',WOPT 
WRITE(11,'(A,A)')' Mean removal: ',RM 
WRITE(11,'(A,A12)')! Input file: ',INFILE 


Start the outer loop, the segment loop **** 
DO 50 NN=1,NSEG 
WRITE(*,*)' Computing averaged PSD from segment #',NN 
initialize spectrum *** 
OO 53° 2=2 Meer) 
P(I)=0. 
TPWR(I)=0. 
PCTPWR(I)=0. 
DBP(I)=0. 
CONTINUE 
DEN=0. 


Read M Points from pulse # L (starting with ISTART) and write to 
ARRAY Wl **x* 
Start the pulse loop *** 
DO 120 L=1,NPLS 
and the data point (time) loop *** 
DO 100 I=1,M 
Tit=I+iSTART=1 
L2=I+1I 
put the real values in the odd numbered elements of Wl. *** 
W1(I2-1) =ENSMBL(II,L) 
the imaginary part (even elements) is zero *** 
W1(I2)=0. 
CONTINUE 


Now remove the mean (of the segment) if it was selected above *** 
IF(RM .EQ. ‘Y* .OR. RM EQ. ‘yj Pen 
SUMRM=0. 
DO 105 I=1,M 
I2=I+I 
SUMRM=SUMRM+W1 (I2-1) 
CONTINUE 
REVMN=SUMRM/M 
DO 106 I=1,M 
T2=I+1 
W1(I2-1)=W1 (12-1) -REVMN 
CONTINUE 
ENDIF 


Apply the window function to the M data points. Operate only on the 
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* real part of Wl since the imag. part 1s Zero. **** 
DO 160 J=1,M 
J2=St+I 
W=WINDOW(J,M, WOPT) 
W1(J2-1)=W1(J2-1) *W 
160 CONTINUE 


wet Fourier transform the windowed data, then sum results into previous 
* segment for later averaging *** 

CALL FOUR1(W1,M,1) 
wae now sum up the real pos. and neg. frequencies from each segment *** 


P(1)=P(1)+W1(1) **2 
DO 17 J=2,MO02 
J2=Jt+I 
P(J)=P(J)+W1(J2<-1) ®**24+W1 (M43-J2) **2 
a7 CONTINUE 
P (MO2+1)=P(MO2+1) +W1 (M+1) **2 
DEN=DEN+SUMW 
ak now, go get a segment from the next pulse and do it again *** 
120 CONTINUE 


C *** correct normalization for windowing *** 
DEN=M*4*DEN 
C ***x normalize the output, ie. compute the average power over all pulses *** 
DO 19 J=1,M02+1 
P(J)=P(J) /DEN 
19 CONTINUE 
* 
Cx** compute the relative PSD in decibels and the running total power *** 
* zero GB is the peak value *** 
PMAX=P(1) 
DO 51 J=1,MO02 
IF(P(J+1) .GT. PMAX) PMAX=P(J+1) 
51 CONTINUE 
SUMPWR=0. 
DO 60 J=1,MO02+1 
TFe UPd) 250.0.) THEN 
WRITE(*,*) 'WARNING: P(J) has a Zero value. P(J) set to -999 
u999260" “at. a J of',J,' cemCOnel ntl nag .. 
DBP (J) =-999999. 
ELSE 
DBP(J)=20. *ALOG10(P(J)/PMAX) 
SUMPWR=SUMPWR+P (J) 
TPWR(J) =SUMPWR 
ENDIF 
60 CONTINUE 
* 
*** compute the percentage of running total power at freq. J relative to the 
* total power summed over all freqs. *** 
DO 70 J=1,M02+1 
PCTPWR(J)=100. *TPWR(J) /TPWR(MO2+1) 
70 CONTINUE 
alialial write out the results for segment # NN *** 
WRITE(11,'(3(A,I6))')' PSD OF TIME SEGMENT',NN,', POINT NO.',IST 
LART, ' THROUGH', ISTART+M-1 
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DO 300 LL=1,MO02+1 
FREQ=LL*DFREQ - DFREQ 
WRITE(11,°(1X,F10.1,F17.3,F10.3,F17.3,F8.2)") FREO pee eee ee 
1LL) , TPWR(LL) , PCTPWR( LL) 
300 CONTINUE 
* 
wae now get ready for the next segment *** 
ISTART=ISTART+IWIDE 
50 CONTINUE 
* 
*** All Done *** 
z*x* If desired, go back and run another operation on the same ensemble *** 
WRITE(*,'(/A\)')' DO YOU WISH TO CONTINUE WITH THIS ENSEMBLE? (Y o 
ircNje 
READ(*,'(A)') CONT 
IF (CONT .EQ. “Y" 2OR. CONT. EO.) 9) )) erHeEN 
CLOSE(11) 
PRINT '(/A\)',' ENTER ANOTHER (NEW) OUTPUT DATA FILE NAME: ! 
READ (*,-' CA) ") OUTFILE 
OPEN (11, FILE=OUTFILE, STATUS='NEW') 
GO TO 29 
ENDIF 
CLOSE (10) 
CLOSE (21) 
END 
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FUNCTION WINDOW(J,MM,WOPT) 


These window functions come from the paper "On the Use of Windows 
for Harmonic Analysis with the Discrete Foutier Transform", by 
Fredric J. Harris, IEEE Proc., Vol. 66,°Ne.1, Jan,1970. 
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INTEGER WOPT 

PI2=8.*ATAN(1.) 

M=MM/2 

FACM=M-0.5 

FACP=1./(M+0.5) 

IF(WOPT .EQ. 1) THEN 
Cc -- PARZEN (TRIANGULAR) WINDOW: 
WINDOW=(1.-ABS(((J-1)-FACM) *FACP) ) 


Cc 
ELSEILF(WOPRT .EO.. 2). sHEN 

c -- RECTANGULAR WINDOW: ie. no window 
WINDOW=1. 

Cc 


ELSEIF(WOPT .EQ. 3) THEN 
e -- WELCH (PARABOLIC) WINDOW: 
WINDOW=(1.-(((J-1)-FACM) *FACP) **2) 


ELSEIF(WOPT .EQ. 4) THEN 


C -- VAN HANN (SQUARED COSINE) WINDOW: 
WINDOW=0.5-0.5*COS(PI2* (FLOAT (J-1) /FLOAT (MM) ) ) 
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BESEIF(WOPT -EQ@: 5) THEN 
C -- HAMMING (RAISED COSINE) WINDOW: 
WINDOW=0.538-0.462%*COS (PI2* (FLOAT (J-1) /FLOAT (MM) ) ) 


c 
ELSEIF(WOPT .EQ. 6) THEN 
Cc -- BLACKMAN-HARRIS 3 TERM MINIMUM WINDOW: 
WINDOW =.42323-.49755*COS (PI2* (FLOAT (J-1) /FLOAT (MM) ))+ 
1 .07922*COS (PI2*2.* (FLOAT (J=-1) /FLOAT (MM) ) ) 
S 
ELSEIF(WOPT .EQ. 7) THEN 
eC -- BLACKMAN-HARRIS 4 TERM MINIMUM WINDOW: 
WINDOW =.35875-.48829*COS (PI2* (FLOAT(J-1) /FLOAT(MM) ))+ 
1 .14128*COS (PI2*2.* (FLOAT (J-1) /FLOAT (MM) ))- 
2 .01168*COS (PI2*3.* (FLOAT (J-1) /FLOAT(MM)) ) 
c 
ELSE 
WRITE(*,*) "INVALID WINDOW SELECTION, TERMINATE...' 
STOP 
ENDIF 
RETURN 
END 


CrexkkkkkkkkkkkakekkrkeeeaekkrkeeeeeeKeeKeKRKRKRKRKRKRKAKRKKKKKKRKRKRKRKRKKKKKRKRKKKKRKKRKRKKRhKK 


SUBROUTINE FOUR1 (DATA,NN,ISIGN) 


FFT subroutine comes from the book "Numerical Recipes", by Press et.al. 
Cambridge Univ. Press, 1986. 
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REAL*8 WR,WI,WPR,WPI,WTEMP, THETA 
DIMENSION DATA(*) 
N=2*NN 
J=1 
DO 11 I=1,N,2 
IF(J.GT.1I) THEN 
TEMPR=DATA (J) 
TEMPI=DATA (J+1) 
DATA (J)=DATA (I) 
DATA (J+1)=DATA(I+1) 
DATA (I)=TEMPR 
DATA (I+1) =TEMPI 
ENDIF 
M=N/2 
il IF ((M.GE.2).AND.(J.GT.M)) THEN 
J=J=-M° 
=M/2 
Go: TO 1 
ENDIF 
J=I+M 
ia CONTINUE 
MMAX=2 
2 IF (N.GT.MMAX) THEN 
ISTEP=2 *MMAX 
THETA=6 . 28318530717959D0/ (ISIGN*MMAX) 
WPR=-2.D0O*DSIN(0.5D0O*THETA) **2 
WPI=DSIN(THETA) 
WR=1.D0 
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wI=0.DO0 
DO 13 M=1,MMAX,2 
DO 12 I=M,N,ISTEP 
J=I+MMAX 
TEMPR=SNGL(WR) *DATA(J) -SNGL(WI) *DATA(J+1) 
TEMPI=SNGL(WR) *DATA(J+1)+SNGL(WI) *DATA(J) 
DATA (J) =DATA (I) -TEMPR 
DATA (J+1) =DATA (I+1) -TEMPI 
DATA (I) =DATA(I)+TEMPR 
DATA (I+1) =DATA (I+1)+TEMPI 
1) CONTINUE 
WTEMP=WR 
=WR*WPR-WI1 *WPI+WR 
WI=WI*WPR+WTEMP*WPI+WI 
13 CONTINUE 
MMAX=ISTEP 
GO To 2 
ENDIF 
RETURN 
END 


Z02 


J. 


USER’S GUIDE FOR PROGRAM CORREX 


1. Input Instructions 


Insure that the executable file for program CORREX 
(CORREX.EXE) and the input data file are in the same 
subdirectory of a disk drive with sufficient space 
remaining to write the output file. Not that the input 
file format for program CORREX depends on whether 
autocorrelation or cross correlation is desired. lige 
autocorrelation is desired, the input file is expected 
to be a sequential file containing samples of a single 
variable (e.g., u component velocity). If cross 
correlation is desired, the input file is expected to be 
a sequential ASCII file containing alternating samples 
of two variables (e.g., u component and v component 
velocity). 


Initiate program CORREX (type "CORREX" at the DOS 
ProOmpe). 


The program will prompt the user for the name of the 
input data file. Enter the full name of the input data 
file including the extension. 


The program will prompt the user for the name of the 
output file which will be created by the program. Enter 
the desired file name. Do not use the name of an 
existing file or the program will terminate. 


The program will prompt the user to select the option for 
autocorrelation (1) or cross correlation (2). Enter the 
desired choice. 


The program will prompt the user to enter the number of 
lags which are desired as a fraction of the total number 
of data points. For example, if the autocorrelation is 
desired on a time series containing 1000 points, and the 
user wishes to compute the autocorrelation for 250 lag 
points, then 0.25 is entered at the prompt. For truly 
random data, the autocorrelation usually falls rapidly 
to near zero such that typical lag values are 0.1 or 0.2 


The program will compute the selected correlation 
function and write the output file. 
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2. Output Description 
The output file consists of two columns. The first 
column contains the lag number. The lag number divided by the 
sample rate yields the value of the lag time, . The second 
column presents the selected correlation function, 
corresponding to the lag numbers. The following pages present 


the source code for program CORREX. 
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$DEBUG 


PROGRAM CORREX 


pe ee eee eee SEES ESE SE ESE SEE SERS REE ERS ESE SESE EES ERS ERE RE RE RE REESE ERS SEER ESSE ES 
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This program computes an estimate of the autocorrelation function 

or cross correlation function of a real discrete (stationary) time series 
using the explicit method. The main program is a driver to set up the 
input for the call to SUBROUTINE CORMAR. Subroutine CORMAR comes from 
"Digital Spectral Analysis with Applications" by S. Lawrence Marple, 
Prentiss-Hall, 1987. The explicit method (see CORMAR) is slower than 

the more typical FFT methods, but 1s easier to understand and progran. 

It also has the advantage of not requiring an integer power 

of 2 number of data points like the FFT versions require. 


For autocorrelation, the input file should be a sequential file containing 
values of velocity (or whatever single parameter you wish to 
autocorrelate). The time variable is implied knowing the sample rate. 


For cross correlation, the input file should be a sequential file 
containing alternating values of the two component velocities at each 
time sample (or whichever parameters you wish). As above, the time is 
implied by the sample rate. 


After reading the data into the proper arrays and calibrating, the mean 

1s computed and subtracted from DATA1l and DATA2 to obtain only the 
fluctuating velocity component u' (and v' for cross correl.). These 

values are stored back in DATA1 and DATA2. Note, this program only 

computes correlation for positive lag values. Autocorrelations are 

even functions, so you don't need the negative lags. For cross correlations 
Simply reverse the input data values so that DATAl1 becomes DATA2 and visa 
versa. This works because Ruv(-tau) = Rvu(tau). 


The key variables are defined below: 


N = Array dimension parameter, maximum number of data points 

NPTS = Number of velocity input values. 

DATA1 = First data array (i.e. u') 

DATA2 = second data array (1.e. v', is set equal to u for autocorrel.) 

LAGG = Maximum dimension of correlation array. 

LAG = The number of lags for which to compute the correlation (less than N) 
RO = The correlation result at lag=0 

R = The resulting correlation array for lags 1 through LAG. 


Programmed by: Donald K. Johnson 
August 1988, NPS Monterey CA. 


tee ee ee eee eRe ERE SEER ESSE ELSES EEE SEP EE SER EP ES EPEC ESE SESE SSS ESS EL ELS SL EL ESS 


PARAMETER (N=4096, LAGG=4096) 
DIMENSION DATA1(N) ,DATA2(N) ,R(LAGG) 
INTEGER NPTS, LAG 

CHARACTER*12 INFILE,OUTFILE 


PRINT '(/A/)‘',' Computes the auto or cross correlation function 
lusing the explicit method. ' 
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PRINT '(/A\)',' ENTER THE INPUT DATA FILE NAME: ' 
READ(*,'(A)') INFILE 

PRINT '(/A\)',' ENTER THE OUTPUT FILE NAME: ' 
READ(*,'(A)') OUTFILE 


OPEN (9,FILE=INFILE,STATUS='OLD') 
OPEN(10, FILE=OUTFILE, STATUS='NEW') 
* 
*** choose the input file format **** 
5 PRINT '(/A)',' Do you want to compute autocorrelation <1> or cro 
lss correlation <2> ?' 
PRINT "(AN)" ENTER CHOICE: == 
READ(*,*) KOPT 


PRINT '(/A\)',' ENTER THE NUMBER OF LAG POINTS (A DECIMAL FRACTI 
1ON OF THE TOTAL NUMBER OF POINTS = usually around 0.1] or eGe 
22) ae 
READ(*,*) FRAC 
® 
kat currently deactivated 
* PRINT '(/A\)',' ENTER THE CALIBRATION COEFFICIENTS A AND B [wher 
* le u = A*(data**B)): ' 
* READ(*,*) A,B 
x 
xt for autocorrelation *** 


IF(KOPT .EQ. 1) THEN 
DO 11 I=1,N 
READ (9,*,END=15,ERR=15) DATA1(I) 


*** convert the hot wire voltage to velocity in ft./sec. *** 


al currently deactivated 
kee DATA1(1I)=A* (DATA1(I) **B) 
* 
DATA2 (I) =DATA1 (I) 
sleet CONTINUE 
* 
eek for cross correlation *** 


ELSELF (KOFI. «EO. 2) ° THEN 
DO 12 I=1,N 
READ(9,*,END=15,ERR=15) DATA1(1) ,DATA2 (TI) 


*** convert the hot wire voltage to velocity in ft./sec. *** 


alia currently deactivated 
* DATA1(I)=A* (DATA1 (I) **B) 
* DATA2 (I) =A* (DATA2 (I) **B) 
x 
LZ CONTINUE 
EESE 
PRINT '(/A/)',* Invalid choice, try again. 
GO TO 5 
ENDIF 
15 CONTINUE 


x 
z*x* compute the maximum lag based on the number of data points *** 
NPTS=I 
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LAG=INT (NPTS*FRAC) 


remove the mean i.e. compute u' (and v') *** 
SUMD1=0. 
SUMD2=0. 

DO 100 I=1,NPTS 
SUMD1=SUMD1+DATA1 (T) 
SUMD2=SUMD2+DATA2 (T) 

CONTINUE 

AVG1=SUMD1/NPTS 

AVG2=SUMD2/NPTS 

DO 200 I=1,NPTS 
DATA1 (I) =DATA1 (I) -AVG1 
DATA2 (I) =DATA2 (I) -AVG2 

CONTINUE 


compute the unbiased correlation estimate **** 


PRINT (7A) -Computing correlation ...* 
CALL CORMAR (NPTS,LAG,0,DATA1,DATA2,RO,R) 


Write the output. Note, this program only computes correlations 
for positive lags. 
IZERO=0 
WRITE(10, '(1X,1I5,3X,F15.3)') IZERO,RO 
DO 13 I=1,LAG 
WRELE (NOs CIN PES (5k, 615.3) °) 1, RC1) 
CONTINUE 
CLOSE(9) 
CLOSE(10) 
END 


KAMARA KKKKKKEKEKREKREKEKERRKARAKEERRERERRRREKKKEERRRRKKRKKKARKERKEK 


SUBROUTINE CORMAR (N,LAG,MODE,X,Y,RO,R) 


Subroutine from S.L. Marple's "Digital Spectral Analysis with 
Applications", Prentiss-Hall, 1987 


Modified for real valued inputs by Donald Johnson, Aug. 1988, Naval 
Postgraduate School, Monterey, CA. 


This program computes either the unbiased or biased real correlation 
estimates between real data sample arrays X and Y. If X=Y, then the 
autocorrelation is computed. 


Input Parameters: 


N ~ Number of data samples in arrays X and Y (integer) 

LAG - Number of correlation lags to compute [{ lags from 0 to LAG 
are computed and stored in RO and R(1) to R(LAG) ] (integer) 

MODE - Set to 0 for unbiased correlation estimates; otherwise, biased 
correlation estimates are computed (integer) 

X - Array of real data samples X(1) through X(N) 

Y - Array of real data samples Y(1) through Y(N) 


Output Parameters: 
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RO - Real correlation estimate for lag 0 
R - Array of real correlation estimates for lags 1 to LAG 
Notes: 


External arrays X,Y must be dimensioned .GE. N and array R must be 
dimensioned .GE. LAG in the calling progran. 


QAAAANAAANANANDA 


REAL X(1),¥(1),R(1),RO,SUM 
INTEGER N,LAG,K,MODE,NK 
DO 30 K=0,LAG 
NK=N~K 
SUM=0. 
DO 10 J=1,NK 
10 SUM=SUM+X (J+K) *Y (J) 
IF (K .NE. 0) GO TO 20 
RO=SUM/ FLOAT (N) 
GO TO 30 
20 IF (MODE .EQ. 0) R(K)=SUM/FLOAT(N-K) 
IF (MODE .NE. 0) R(K)=SUM/FLOAT(N) 
30 CONTINUE 
RETURN 
END 
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USER’S GUIDE FOR PROGRAM CORFFT 


1. Input Instructions 


Insure that the executable file for program CORFFT 
(CORFFT.EXE) and the input data file are in the same 
subdirectory of a disk drive with sufficient space 
remaining to write the output file. Not that the input 
file format for program CORFFT depends on whether 
autocorrelation or cross correlation is desired. If 
autocorrelation is desired, the input file is expected 
to be a sequential file containing samples of a single 
variable (e.g., wu component velocity). If cross 
correlation is desired, the input file is expected to be 
a sequential ASCII file containing alternating samples 
of two variables (e.g., u component and v component 
velocity). 


Initiate program CORFFT (type "CORFFT' at the DOS 
prompt). 


The program will prompt the user for the name of the 
input data file. Enter the full name of the input data 
file including the extension. 


The program will prompt the user for the name of the 
output file which will be created by the program. Enter 
the desired file name. Do not use the name of an 
existing file of the program will terminate. 


The program will prompt the user to select the option for 
autocorrealtion (1) or cross correlation (2). Enter the 
number of the desired choice. 


The program will prompt the user to enter the number of 
lags which are desired as a fraction of the total number 
of data points. For example, if the autocorrelation is 
desired on a time series containing 1000 points, and the 
user wishes to compute the autocorrelation for 250 lag 
points, then 0.25 is entered at the prompt. For truly 
random data, the autocorrelation usually falls rapidly 
to near zero such that typical lag values are 0.1 or 0.2. 


The program will compute the selected correlation 
function and write the output file. 
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2. Output Description 
The output file consists of two columns. The first 
column contains the lag number. The lag number divided by 
the sample rate yields the value of the lag time, . The 
second column presents the selected correlation function, 
corresponding to the lag numbers. The following pages present 


the source code for program CORFFT. 
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$ DEBUG 


PROGRAM CORFFT 


KAKA ARARAAKAARAARAEEAAKARAEAREREREAEAAKKARREAAKRREAAKAKEREAEEKEAEKKKRKRKAKRAKAEE 
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This program computes an estimate of the autocorrelation function 

or cross correlation function of a real discrete (stationary) time 
series using the Fourier transform method. The main program is a driver 
to set up the input for the call to SUBROUTINE CORREL. Subroutine CORREL, 
TWOFFT, REALFT AND FOUR] come from "Numerical Recipes" by Press et. al., 
Cambridge, 1986. 


For autocorrelation, the input file should be a sequential file containing 
values of velocity (or whatever parameter you wish to autocorrelate). 
The time variable is implied knowing the sample rate. 


For cross correlation, the input file should be a sequential file 
containing alternating values of the two component velocities at each 
time sample (or whichever parameters you wish). AS above, the time is 
implied by the sample rate. 


After reading the data into the proper arrays and calibrating, the mean 
is computed and subtracted from DATAl1 and DATA2 to obtain only the 
fluctuating velocity component u' (and v' for cross correl.). These 
values are stored back in DATA1 and DATA2. 


The Fourier transform method of computing correlation results ina 
circular correlation. Two convert this to linear correlation, the elements 
of array ANS are divided by (NPTS-I) where I is the current lag value. I 
varies from 0 to LAG. This program is only outputting correlation values 
at positive lags; ANS(1) for zero lag to ANS(N/2) for lag=LAG. 
Autocorrelations of real time series are even functions so yopu don't 
need the negative lags. For cross correlations you may want to see the 
negative lags, so just modify the output WRITE statements or reverse the 
input data vectors. Reversing the input time series will give you the 
negative lag values because Ruv(-tau) = Rvu(tau). Increasingly 

negative lags are stored in wraparound mode from ANS(N) for lag -1 

to ANS(N-LAG+1) for -LAG. 


The key variables are defined below: 


N = Array dimension parameter, maximum number of data points 

NPTS = Number of velocity input values. 

= First data array (l.e. u') 

DATA2 = second data array (i.e. v', is set equal to u for autocorrel.) 

INT2 = The size of DATAl1 and DATA2 after zero padding to the smallest 
integer power of 2 greater than NPTS+LAG. (The size of the FFT.) 

The number of lags for which to compute the correlation (less than N) 

The circular correlation result for lags 0 through LAG. 


E 
@ 


Programmed by: Donald K. Johnson 
August 1988, NPS Monterey CA. 
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PARAMETER (N=4096, N2=2*N) 
DIMENSION DATA1(N) , DATA2(N) ,ANS(N2) 
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INTEGER NPTS, LAG 
CHARACTER*12 INFILE, OUTFILE 


PRINT '(/A/)',' Computes the auto or cross correlation function 
lusing Fourier transform methods.' 


PRINT '(/A\)',' ENTER THE INPUT DATA FILE NAME: ' 
READ(*,'(A)') INFILE 
PRINT '(/A\)',' ENTER THE OUTPUT FILE NAME: ' 


READ(*,'(A)') OUTFILE 


OPEN (9, FILE=INFILE,STATUS='OLD') 
OPEN(10, FILE=OUTFILE, STATUS= NEW") 


choose the input file format **** 


PRINT '(/A)',' Do you want to compute autocorrelation <1> or cro 
lss correlation <2> ?! 

PRINE “GAN ENTER GHOTCE:. = 

READ (*,*) KOPT 

PRINT '(/A\)',' ENTER THE NUMBER OF LAG POINTS (A DECIMAL FRACTI 
10N OF THE TOTAL NUMBER OF POINTS - usually around 0. 2 o7MG7 


Ze 
READ(*,*) FRAC 


currently deactivated 

PRINT '(/A\)',' ENTER THE CALIBRATION COEFFICIENTS A AND B [wher 
le u = A*(data**B)j): ' 

READ(*,*) A,B 


for autocorrelation *** 
IF(KOPT .EQ. 1) THEN 
DO 11 I=1,N 
READ (9,*,END=15,ERR=15) DATA1(TI) 


convert the hot wire voltage to velocity in ft./sec. *** 
currently deactivated 
DATA1 (I) =A* (DATA1 (I) **B) 


DATA2 (I) =DATA1 (T) 
CONTINUE 


for cross correlation *** 
ELSEIF (KOPT ».EO.. 2) CHER 
DO 12 I=1,N 
READ (9, *,END=15,ERR=15) DATA1(I) , DATA2(TI) 


convert the hot wire voltage to velocity in ft./sec. *** 
currently deactivated 

DATA1(I)=A*(DATA1 (I) **B) 

DATA2 (I) =A* (DATA2 (I) **B) 


CONTINUE 


ELSE 
PRINT '(/A/)',' Invalid choice, try again.' 
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GO TO 5 
ENDIF 
CONTINUE 


compute the maximum lag based on the number of data points *** 
NPTS=I 
LAG=INT (NPTS * FRAC) 


zero pad the data arrays with at least LAG zeros such that the data 
array is an integer power of two. Max array size is currently 4096. 
DO 25 J=1,12 
INT2=2**J 
find the smallest power of 2 array size bigger than NPTS+LAG *** 
IF(INT2 .GT. NPTS+LAG) GO TO 27 
CONTINUE 
CONTINUE 
zero pad from NPTS to INT2 
DO 30 I=NPTS+1,INT2 


DATA1 (I) =0. 
DATA2(I)=0. 
CONTINUE 
remove the mean 1.e. compute u' (and v') *** 
SUMD1=0. 
SUMD2=0. 


DO 100 I=1,NPTS 
SUMD1=SUMD1+DATA1 (I) 
SUMD2=SUMD2+DATA2 (I) 

CONTINUE 

AVG1=SUMD1/NPTS 

AVG2=SUMD2/NPTS 

DO 200 I=1,NPTS 
DATA1(1I)=DATA1(I) -AVG1 
DATA2 (I) =DATA2 (I) -AVG2 

CONTINUE 


compute the unbiased correlation estimate **** 
PRint = (/A\)",° Computing correlation with FFT method...’ 
CALL CORREL (DATA1,DATA2,INT2,ANS) 


Write the output. Note: ANS is divided by NPTS-I to convert the 
Circular correlation to linear correlation. Also, only the positive 
lags are output. To get the increasingly negative lags, run a DO 
loop backwards from NPTS (lag of -1) to NPTS-LAG+1 (lag of -LAG). 


DO 13 I=0,LAG~-1 
WRITE (10, '(1X,I5,3X,F15.3)') I,ANS(I+1)/(NPTS-I) 

CONTINUE 

CLOSE (9) 

CLOSE(10) 

END 
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SUBROUTINE CORREL(DATA1,DATA2,N,ANS) 


* From "Numerical Recipes" by Press et.al., Cambridge, 1986. 
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Computes the correlation of two real data sets DATA1 and DATA2, each of 
length N (including zero padding). N must be an integer power of 2. The 
answer is returned as the first N points in ANS stored in wraparound 
order. (Increasingly negative lags are in ANS(N) down to ANS(N/2+1), 


be at least length 2*N in the calling program since it is also used as 


at positive lags. 


+ + + ee OF Oe Oe OF 


* 
* 


maximum FFT size *** 

PARAMETER (NMAX=8192) 

DIMENSION DATA1(N) ,DATA2(N) 

COMPLEX FFT (NMAX) , ANS (N) 

Lokal Fourier transform both data vectors at once *** 
CALL TWOFFT (DATA1,DATA2,FFT,ANS,N) 

wae normalize for inverse FFT *** 

NO2=FLOAT(N)/2.0 

DO 11 I=1,N/2+1 


& & multiply to find the FFT of their correlation *** 
ANS (I) =FFT (I) *CONJG (ANS (I) ) /NO2 
all CONTINUE 
ANS (1) =CMPLX (REAL(ANS(1)),REAL(ANS(N/2+1) )) 
ee inverse transform to obtain circular correlation *** 
CALL REALFT(ANS,N/2,-1) 
RETURN 
END 


CHR ERA KREEERAAREKE KERR EEKEKREKKAKKKKKKKKKKKKKK 
SUBROUTINE TWOFFT(DATA1,DATA2,FFT1,FFT2,N) 
* 
* From "Numerical Recipes" by Press et.al., Cambridge, 1986. 
* 
DIMENSION DATA1(N) , DATA2(N) 
COMPLEX PEIL(UN) FFI2Z (NB HZ, Cl C2 
C1=CMPLX(0.5,0.0) 
C2=CMPLX(0.0,-0.5) 
DO 11 J=1,N 
FFT1 (J) =CMPLX (DATA1 (J) , DATA2(J) ) 
i bll CONTINUE 
CALL FOUR] (FFT1,N, 1) 
FFT2 (1) =CMPLX (AIMAG(FFT1(1)),0.0) 
FFT1(1)=CMPLX (REAL(FFT1(1)),0.0) 
N2=N+2 
DO 12 J=@2,N/271 
H1=C1*(FFT1(J)+CONJG(FFT1(N2-J) )) 
H2=C2* (FFT1(J) -CONJG(FFT1(N2-J))) 
FFT1(J)=H1 
FFT1(N2-J) =CONJG(H1) 
FFT2 (J) =H2 
FFT2 (N2-J) =CONJG (H2) 
12 CONTINUE 
RETURN 
END 


CRAMER HKERAEARAKEREKEKEKREKEEEAKEEAAEEKEKKKEEKEEKEEEKKKEKK 


SUBROUTINE REALFT (DATA,N, ISIGN) 
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while increasingly positive lags are in ANS(1) up to ANS(N/2) ). ANS must 


workspace. Sign convention: if DATA1 lags DATA2 then ANS will show a peak 
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From "Numerical Recipes" by Press et.al., 


REAL*8 WR,WI,WPR,WPI,WTEMP, THETA 
DIMENSION DATA(*) 
THETA=6. 28318530717959D0/2.0D0/DBLE(N) 
C1=0.5 
IF (ISIGN.EQ.1) THEN 
C2=-0.5 
CALL FOUR1(DATA,N,+1) 
ELSE 
C2=0.5 
THETA=-THETA 
ENDIF 
WPR=-2.0D0*DSIN(0.5DO*THETA) **2 
WPI=DSIN (THETA) 
WR=1.0D0+WPR 
WI=WPI 
N2P3=2*N+3 
DO 11 I=2,N/2+1 
I1=2*I-1 
I2=I1+1 
I3=N2P3-I2 
I4=I3+1 
WRS=SNGL(WR) 
WIS=SNGL(WI) 
H1R=C1*(DATA(I1)+DATA(I3)) 
H1I=C1* (DATA(I2) -DATA(I4)) 
H2R=-C2* (DATA(I2)+DATA(I4)) 
H2I=C2* (DATA (I1)-DATA(I3)) 
DATA (I1)=H1R+WRS *H2R-WIS*H2I 
DATA (12) =H1I+WRS*H2I+WIS*H2R 
DATA (13) =H1R-WRS *H2R+WIS*H2I 
DATA (14) =-H1I+WRS*H2I+WIS*H2R 
WIEMP=WR 
WR=WR*WPR-WI*WPI+WR 
WI=WI*WPR+WTEMP*WPI+WI 
CONTINUE 
IF (ISIGN.EQ.1) THEN 
H1R=DATA(1) 
DATA (1) =H1R+DATA(2) 
DATA (2) =H1R-DATA (2) 
ELSE 
H1R=DATA(1) 
DATA(1)=C1*(H1R+DATA(2)) 
DATA (2) =C1* (H1R-DATA(2)) 
CALL FOUR1(DATA,N,-1) 
ENDIF 
RETURN 
END 


Cambridge, 1986. 
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SUBROUTINE FOUR1(DATA,NN,ISIGN) 
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From "Numerical Recipes" by Press et.al., Cambridge, 1986. 


ng 


iz 
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REAL*8 WR,WI,WPR,WPI,WTEMP, THETA 
DIMENSION DATA(*) 
N=2*NN 
J=1 
DO 11 I=1,N,2 
IF(J.GT.I) THEN 
TEMPR=DATA (J) 
TEMPI=DATA (J+1) 
DATA(J) =DATA(T) 
DATA (J+1) =DATA (I+1) 
DATA(I)=TEMPR 
DATA (I+1) =TEMPI 
ENDIF 
M=N/2 
IF ((M.GE.2).AND.(J.GT.M)) THEN 
J=JI-M 
M=M/2 
GO TO 1 
ENDIF 
J=I+M 
CONTINUE 
MMA X=2 
IF (N.GT.MMAX) THEN 
ISTEP=2 *MMAX 
THETA=6. 28318530717959D0/ (ISIGN*MMAX) 
WPR=-2.DO*DSIN(0.5DO*THETA) **2 
WPI=DSIN (THETA) 
WR=1.D0 
WI=0.DO0 
DO 13 M=1,MMAX,2 
DO 12 I=M,N,ISTEP 
J=I+MMAX 
TEMPR=SNGL(WR) *DATA(J) -SNGL(WI) *DATA(J+1) 
TEMPI=SNGL(WR) *DATA(J+1)+SNGL(WI) *DATA (J) 
DATA (J)=DATA(I) -TEMPR 
DATA (J+1) =DATA(I+1) -TEMPI 
DATA (I) =DATA(I)+TEMPR 
DATA (I+1)=DATA(I+1)+TEMPI 
CONTINUE 
WTEMP=WR 
WR=WR*WPR-WI*WPI+WR 
WI=WI *WPR+WTEMP*WPI+WI 
CONTINUE 
MMAX=ISTEP 
GO TO 2 
ENDIF 
RETURN 
END 
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L. SOURCE CODE FOR PROGRAM MEMPSD 


The following pages present the source code for program 
MEMPSD. This program is included for the sake of completeness 


but without documentation. 
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PROGRAM MEMPSD 


Driver for routine MEMCOF and EVLMEM - Maximum Entropy Method 


PARAMETER (N=16384 ,M=100,NFDT=512) 
DIMENSION DATA(N) ,COF(M) ,WK1(N) ,WK2(N) ,WKM(M) 
identify input and output files 
OPEN (5,FILE='RECORD1.DAT',STATUS='OLD') 
OPEN (7, FILE='PSDOUT1.DAT' ,STATUS='NEW') 
read input file 
READ(5,*) (DATA(I),PROP,I=1,N) 
CLOSE(S) 
compute coefficients of spectral estimate 
CALL MEMCOF(DATA,N,M, PM, COF,WK1,WK2,WKM) 
write output 
WRITE(7,*) ‘Power spectrum estimate of data in RECORD1.DAT' 
WRITE(7,'(1X,T6,A,T25,A,T40,A)°) *f*delta”. freq") po.en. 
SMPLRT=30000. 
DO 11 J=0,NFDT 
EDT =0.5*1 /NEDI 
FREQ=FDT*SMPLRT 
function EVLMEM evaluates the spectral density function using the 
coefficients computed with MEMCOF. 
WRITE(7,' (3 (2X, F14.4))*) FDT, FREQ, EVUMEM (PDT Cor Mer 
CONTINUE 
END 
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SUBROUTINE MEMCOF(DATA,N,M, PM, COF,WK1, WK2, WKM) 


From the book "Numerical Recipes" by Press, et al., Cambridge University 


press, 1986. See this book for documentation. 


DIMENSION DATA(N) ,COF(M) ,WK1(N) ,WK2(N) , WKM(M) 
P=0. 
BOUL J=1,N 

P=P+DATA(J) **2 

CONTINUE 
PM=P/N 
WK1 (1) =DATA(1) 
WK2 (N-1)=DATA (N) 
DO 12 J=2,N-1 
WK1 (J) =DATA(J) 
WK2 (J-1) =DATA(J) 
CONTINUE 
DO 17 K=1,M 

PNEUM=0. 

DENOM=0. 

DO 13 J=1,N-K 
PNEUM=PNEUM+WK1 (J) *WK2 (J) 
DENOM=DENOM+WK1 (J) **2+WK2 (J) **2 

CONTINUE 

COF(K)=2.*PNEUM/DENOM 

PM=PM* (1.-COF(K) **2) 

IF(K.NE.1) THEN 
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14 


15 


16 
17 


DO 14 I=1,K-1 
COF (I) =WKM (I) -COF(K) *WKM(K-T) 
CONTINUE 
ENDIF 
IF(K.EQ.M) RETURN 
DO 15 I=1,K 
WKM (I) =COF(T) 
CONTINUE 
DO 16 J=1,N-K-1 
WK1 (J) =WK1 (J) -WKM(K) *WK2 (J) 
WK2 (J) =WK2 (J+1) -WKM(K) *WK1 (J+1) 
CONTINUE 
CONTINUE 
PAUSE ‘never get here' 
END 
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HL 


FUNCTIONSEVEMEM(FDT, COF,M, PM) 


From the book "Numerical Recipes" by Press, et al., Cambridge University 
press, 1986. See this book for documentation. 


DIMENSION COF(M) 
REAL*8 WR,W1,WPR,WPI,WTEMP, THETA 
THETA=6.28318530717959D0*FDT 
WPR=DCOS (THETA) 
WPI=DSIN (THETA) 
WR=1.D0 
WI=0. DO 
SUMR=1. 
SUMI=0. 
DO 11 I=1,M 
WIEMP=WR 
WR=WR*WPR-WI *WPI 
WI=WI *WPR+WTEMP*WPI 
SUMR=SUMR-COF (1) *SNGL(WR) 
SUMI=SUMI-COF (I) *SNGL(WI) 
CONTINUE 
EVLMEM=PM/ (SUMR**2+SUMI**2) 
RETURN 
END 
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